Matlab最小二乘法曲线拟合的应用实例
曲线拟合的线性最小二乘法及其MATLAB程序
1 曲线拟合的线性最小二乘法及其MATLA程序例7.2.1 给出一组数据点(人,yj列入表7^2中,试用线性最小二乘法求拟合曲线,并用(7.2), (7.3 )和(7.4)式估计其误差,作出拟合曲线表7 - 2例7.2.1的一组数据(x^yj解(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.5068.04];plot(x,y, 'r*'),legend('实验数据(xi,yi)' )xlabel( 'x' ), ylabel( 'y'),title('例7.2.1的数据点(xi,yi) 的散点图’)运行后屏幕显示数据的散点图(略)(3)编写下列MATLAB程序计算f(x)在(x i, y i)处的函数值,即输入程序>> 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.A3+ a2.*x.A2+ 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*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*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4];fy=fi-y; fy2=fy.A2; J=sum(fy.A2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)A2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)A2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)A2+(a4+91/10)A2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)A2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)A2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)A2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)A2Q I为求a1,a2,a3,a4使J达到最小,只需利用极值的必要条件——-0 (k =1,2,3,4),得到关于a1,a2,a3,a4的线性方程组,这可以由下面的MATLAB程序完成,即输入程序>> syms al a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-4913/1000*a1 +289/100*a2-17/10*a3+a4...+171/2)A2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)A2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)A2+(a4+91/10)A2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)A2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)A2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)A2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)A2;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的偏导数如下Ja1仁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/625 Ja31 = 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解线性方程组Jan =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];C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*xA3-7988544102557579/562949953421312*xA2 +1804307491277693/281474976710656*x-4648521160813215/562949953421312故所求的拟合曲线为3 2f(x) =5.0911 x -14.1905 x 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=len gth(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); fy2=fy.A2; Ew=max(fy),E仁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 j,%)与拟合函数f的最大误差Ev,平均误差E和均方根误差E 及其数据点(X j,yj和拟合曲线y=f(x)的图形(略).Ew = E1 = E2 =3.105 4 0.903 4 1.240 9a*exp(71/10*b), a*exp(9/2*b), a*exp(13/5*b), a*exp(3/2*b), a,b 的线性方程组,7.3函数m(x)的选取及其MATLA 程序例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.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 i , y i )处的函数值,即输入程序>> 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), 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(34/5*b), a*exp(51/10*b), a*exp(18/5*b), a*exp(17/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(27/10*b), a*exp(18/5*b)]; fy=fi-y; fy2=fy.A 2; J=sum(fy.A 2)运行后屏幕显示误差平方和如下J =(a*exp(17/2*b)-22963/50)A2+(a*exp(87/10*b)-5281/100)A2+( a*exp(71/10*b)-19827/100)A2+(a*exp(34/5*b)-828/5)A2+(a*exp(51/1 0*b)-5917/100)A2+(a*exp(9/2*b)-2083/50)A2+(a*exp(18/5*b)-648/25 F2+(a*exp(17/5*b)-2237/100)A2+(a*exp(13/5*b)-1347/100)A2+(a*ex p(5/2*b)-1287/100)A2+(a*exp(21/10*b)-1187/100)A2+(a*exp(3/2*b)- 669/100F2+(a*exp(27/10*b)-1487/100)A2+(a*exp(18/5*b)-1211/50)A 2为求a,b 使J 达到最小,只需利用极值的必要条件,得到关于 这可以由下面的 MATLAB 程序完成,即输入程序>> syms a bJ=(a*exp(17/2*b)-22963/50)A2+(a*exp(87/10*b)-5281/100)A2+(a*exp(71/10*b)-19827/100)A2+(a*exp(34/5*b)-828/5)A2+(a*exp(51 /10*b)-5917/100)A2+(a*exp(9/2*b)-2083/50)A2+(a*exp(18/5*b)-648/ 25)A2+(a*exp(17/5*b)-2237/100)A2+(a*exp(13/5*b)-1347/100F2+(a*exp(5/2*b)-1287/100)A2+(a*exp(21/10*b)-1187/100)A2+(a*exp(3/2*b )-669/100)A2+(a*exp(27/10*b )-1487/100F2+(a*exp(18/5*b)-1211/50 )A 2;Ja=diff(J,a); Jb=diff(J,b);Ja仁simple(Ja), Jb仁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)A2+8500*a*exp(17/2*b)A2+6800*a*exp(34/5*b)A2-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)A2+5100*a*exp(51/10*b)A2+4500*a*exp(9/2*b)A2+7200*a*exp(18/5*b)A2+3400*a*exp(17/5*b)A2+2600*a*exp(13/5*b)A2+2500*a*exp(5/2*b)A2+1500*a*exp(3/2*b)A2+2700*a*exp( 27/10*b)A2+8700*a*exp(87 /10*b)A2)用解二元非线性方程组的牛顿法的MATLAB^序求解线性方程组J ai =0,J bi =0,得a = b=2.811 0 0.581 6故所求的拟合曲线(7.13)为f(X)= 2.811 0e^.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.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=le ngth(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.A2;Ew=max(fy),E仁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) 的图形')运行后屏幕显示数据(x i , y i)与拟合函数f的最大误差E w = 390.141 5 ,平均误差E1=36.942 2 和均方根误差E2 = 106.031 7 及其数据点(X i , y i )和拟合曲线y=f( x)的图形(略).7.4 多项式拟合及其MATLA程序例7.4.1给出一组数据点(X i, y i )列入表7 43中,试用线性最小二乘法求拟合曲线, 并用(7.2), (7.3 )和(7.4)式估计其误差,作出拟合曲线.表7- 3例7.4.1的一组数据&』)解(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.1219.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.3721 9.1382故拟合多项式为2f(x) = 2.830 2x -7.372 1x 9.138 2 .(4)编写下面的MATLA龍序估计其误差,并做出拟合曲线和数据的图形.输入程序>>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=le ngth(xi); f= 2.8302 .*xi.A2 -7.3721 .*xi+ 9.1382x=-2.9:0.001:3.6;F= 2.8302 .*x.A2 -7.3721 .*x+8.79;fy=abs(f-y); fy2=fy.A2; 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) 的图形’)运行后屏幕显示数据(X i,yj与拟合函数f的最大误差曰,平均误差E1和均方根误差E2 及其数据点(xj i)和拟合曲线y=f(x)的图形(略).Ew = E1 = E2 =0.745 7, 0.389 2, 0.436 37.5 拟合曲线的线性变换及其MATLA程序例7.5.1 给出一组实验数据点(xpyj的横坐标向量为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=(359.26 165.60 59.1741.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^序画出散点图.在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' ), ylabel( 'y'),title('例7.5.1的数据点(xi,yi) 的散点图’)运行后屏幕显示数据的散点图(略)(2)根据数据散点图,取拟合曲线为y = a e bx(a 0,b = 0) , (7.19)其中a, b是待定系数.令Y = ln y, A = ln a, 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.26 165.60 59.17 41.66 25.92 22.37 13.4712.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.A2;Ew=max(fy), E仁sum(fy)/n,E2=sqrt((sum(fy2))/n)运行后屏幕显示y =a e bx的系数b =0.624 1 , a =2.703 9, 数据(x i, y i)与拟合函数f的最大误差Ew =67.641 9 ,平均误差E1=8.677 6 和均方根误差巳=20.711 3 及其数据点(x i,y i)和拟合曲线f (x) =2.703 9e06241x的图形(略).7.6 函数逼近及其MATLA程序最佳均方逼近的MATLABfc程序fun ction [yy1,a,WE]=zjjfbj(f,X,Y,xx)m=size(f); n=le ngth(X);m=m(1);b=zeros(m,m); c=zeros(m,1); if n ~=le ngth(Y) error( 'X和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));endendc(j)=0;for i=1: nc(j)=c(j)+feval(f(j,:),X(i))*Y(i);endenda=b\c;WE=0;for i=1: nff=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:ml=[];for j=1:le ngth(xx)l=[l,feval(f(i,:),xx(j))];endyy=[yy l'];endyy=yy*a; yy1=yy'; a=a';WE;例761 对数据X和Y,用函数y =1, y = x, y =X2进行逼近,用所得到的逼近函数计算在X=6.5处的函数值,并估计误差•其中X=(1 3 4 5 6 7 8 9); Y=(-11 -13 -11 -7 -1 7 17 29).解在MATLA工作窗口输入程序>> X=[ 1 3 4 5 6 7 8 9]; Y=[-11 -13 -11 -7 -1 7 1729];f=['funO';'fun1';'fun2']; [yy,a,WE]=zjjfbj(f,X,Y,6.5)运行后屏幕显示如下yy =2.75000000000003a =-7.00000000000010 -4.99999999999995 1.00000000000000WE =7.172323350269439e-027例7.6.2 对数据X 和Y,用函数y=1, y = x, y = x2,y = cosx,y = e x,y = sin x 进行逼近,其中X =( 0 0.50 1.00 1.50 2.00 2.50 3.00 ),Y=( 0 0.4794 0.8415 0.98150.9126 0.5985 0.1645 ).解在MATLA工作窗口输入程序>> 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=['fu n0';'fu n1';'fu n2';'fu n3';'fu n4';'fu n5'];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.9236Colu mns 8 through 140.9771 0.9926 0.9691 0.9069 0.8080 0.67660.5191Colu mns 15 through 160.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*x A2+0.0765*exp(x) -0.4598*cos(x) +0.5653*s in(x)7.7 三角多项式逼近及其MATLA程序计算三角多项式的MATLA吐程序fun ction [A,B,Y1,Rm]=sanjiao(X,Y,X1,m)n= len gth(X)-1;max1=fix(( n-1)/2);if m > max1m=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)=s in (i*X)*Y'; A(i+1)=cos(i*X)*Y';endA=2*A/n; B=2*B/n; A(1)=A(1)/2;Y 1=A(1);for k=1:m丫仁 Y1+A(k+1)*cos(k*X1)+ B(k+1)*s in (k*X1);Tm=A(1)+A(k+1).*cos(k*X)+ B(k+1).*si n( k*X); k=k+1;endY;Tm; Rm=(sum(Y-Tm).A2)/n;2 i下例7.7.1 根据[-二,二]上的n -13, 60, 350个等距横坐标点X j二-—旦一nX(i =01,2「, n)和函数f(x)=2sin .3(1 )求f(x)的6阶三角多项式逼近,计算均方误差;(2)将这三个三角多项式分别与 f (x)的傅里叶级数一、18.3二“I n f (x) (_1) 2si nnxn y 9n _1的前6项进行比较;(3)利用三角多项式分别计算X i= -2, 2.5的值;(4)在同一坐标系中,画出函数 f (x),n =13,60, 350的三角多项式和数据点的图形•解(1 )输入程序>> X仁-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*si n(X2/3);[A2,B2,Y12,Rm2]=sanjiao(X2,Y2,X1i,6)X3=-pi:2*pi/350:pi;Y3=2*si n(X3/3);[A3,B3,Y13,Rm3]=sanjiao(X3,Y3,X1i,6)X1i=[-2,2.5];Y1=2*si n(X1i/3)for n=1:6bi=(-1)A( n+1)*18*sqrt(3)* n/(pi*(9* n A2-1))end(2)画图,输入程序>>X1=-pi:2*pi/13:pi;Y1=2*si n(X1/3);Xi=-pi:0.001:pi; f=2*si n(Xi/3); [A1,B1,Y1i,R1m]=sanjiao(X1,Y1,Xi,6);X2=-pi:2*pi/60:pi;Y2=2*si n(X2/3); X3=-pi:2*pi/350:pi;Y3=2*s in (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')lege nd('数据点(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 随机数据点上的二元拟合及其MATLA程序例7.8.1设节点(X,Y,Z )中的X和Y分别是在区间[-3, 3]和[-2.5, 3.5]上的50个2 2随机数,Z是函数Z=7-3 x3e-x -y在(X,Y )的值,拟合点(X I ,Y丨)中的X I =-3:023, Y I=-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基三次内插法和MATLAB 4网格化坐标方法计算在(X I,Y I)处的值,作出它们的图形,并与被拟和曲面进行比较.解 (1)最近邻内插法.输入程序>> x=ra nd(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.A3 .* exp(-X.A2 - 丫.人2); %在每个随机点(X,Y )处计算Z的值.X1=-3:0.2:3;%将坐标(XI,YI )网格化.'nearest' ) %计算在每个插值点(XI,YI )xlabel(title( 和节点的图形’)%作二元拟合图形. 'y' ), zlabel( 'z'),mesh(XI,YI, ZI)'x' ), ylabel(用最近邻内插法拟合函数z =7-3 x A3 exp(-x A2 - y A2) 的曲面hold onplot3(X,Y,Z, 'bo') hold of运行后屏幕显示用最近邻内插法拟合函数插值乙(略).(2)三角基线性内插法.输入程序%在当前图形上添加新图形•%用兰色小圆圈画出每个节点(X,Y,Z).%结束在当前图形上添加新图形•2 2Z=7-3 x3e-x -y在两组不同节点处的曲面及其%将坐标(XI,YI )网格化•'lin ear' ) %计算在每个插值点(XI,YI )%作二元拟合图形'y' ), zlabel( 'z'),title(曲面和节点的图形用三角基线性内插法拟合函数)z =7-3 xA3 exp(-xA2 - 丫人hold onplot3(X,Y,Z, 'bo') hold of运行后屏幕显示用三角基线性内插法拟合函数和节点的图形及其插值乙(略)•(3)三角基三次内插法•输入程序%在当前图形上添加新图形•%用兰色小圆圈画出每个节点(X,Y,Z). %结束在当前图形上添加新图形•2 2Z=7-3 x3e-x-y在两组不同节点处的曲面%利用y生成上的随机变量•%在每个随机点(X,Y )处计算Z%将坐标(XI,YI )网格化•'cubic' ) %计算在每个插值点(XI,YI )%作二元拟合图形•),zlabel( 'z'),丫1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1);ZI=griddata(X,Y,Z,XI,YI,处的插值ZI.%legend('拟合曲面','节点(xi,yi,zi)' )>> 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.A3 .* exp(-X.A2 - 丫人2); %在每个随机点(X,Y )处计算Z 的值•X1=-3:0.2:3;Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1);ZI=griddata(X,Y,Z,XI,YI, 处的插值ZI.mesh(XI,YI, ZI)xlabel( 'x' ), ylabel(%legend('拟合曲面','节点(xi,yi,zi)' )>> 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;Z=7-3* X43 •* exp(-X42 - 丫。
曲线拟合的最小二乘法实验
Lab04.曲线拟合的最小二乘法实验【实验目的和要求】1.让学生体验曲线拟合的最小二乘法,加深对曲线拟合的最小二乘法的理解;2.掌握函数ployfit和函数lsqcurvefit功能和使用方法,分别用这两个函数进行多项式拟合和非多项式拟合。
【实验内容】1.在Matlab命令窗口,用help命令查询函数polyfit和函数lsqcurvefit 功能和使用方法。
2.用多项式y=x3-6x2+5x-3,产生一组数据(xi,yi)(i=1,2,…,n),再在yi上添加随机干扰(可用rand产生(0,1)均匀分布随机数,或用randn产生N(0,1)均匀分布随机数),然后对xi和添加了随机干扰的yi用Matlab提供的函数ployfit用3次多项式拟合,将结果与原系数比较。
再作2或4次多项式拟合,分析所得结果。
3.用电压V=10伏的电池给电容器充电,电容器上t时刻的电压为,其中V0是电容器的初始电压,τ是充电常数。
对于下面的一组t,v数据,用Matlab提供的函数lsqcurvefit确定V0和τ。
t(秒) 0.5 1 2 3 4 5 7 9v(伏) 6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63 【实验仪器与软件】1.CPU主频在1GHz以上,内存在128Mb以上的PC;2.Matlab 6.0及以上版本。
实验讲评:实验成绩:评阅教师:200 年月日问题及算法分析:1、利用help命令,在MATLAB中查找polyfit和lsqcurvefit函数的用法。
2、在一组数据(xi,yi)(i=1,2,…,n)上,对yi上添加随机干扰,运用多项式拟合函数,对数据进行拟合(分别用2次,3次,4次拟合),分析拟合的效果。
3、根据t和V的关系画散点图,再根据给定的函数运用最小二乘拟合函数,确定其相应参数。
第一题:(1)>> help polyfitPOLYFIT Fit polynomial to data.P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) ofdegree N that fits the data Y best in a least-squares sense. P is arow vector of length N+1 containing the polynomial coefficients indescending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).[P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and astructure S for use with POLYVAL to obtain error estimates forpredictions. S contains fields for the triangular factor (R) from a QRdecomposition of the Vandermonde matrix of X, the degrees of freedom(df), and the norm of the residuals (normr). If the data Y are random,an estimate of the covariance matrix of P is(Rinv*Rinv')*normr^2/df,where Rinv is the inverse of R.[P,S,MU] = POLYFIT(X,Y,N) finds the coefficients of a polynomial inXHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). Thiscentering and scaling transformation improves the numerical propertiesof both the polynomial and the fitting algorithm.Warning messages result if N is >= length(X), if X has repeated, ornearly repeated, points, or if X might need centering and scaling.Class support for inputs X,Y:float: double, singleSee also poly, polyval, roots.Reference page in Help browserdoc polyfit>>(2)>> help lsqcurvefitLSQCURVEFIT solves non-linear least squares problems.LSQCURVEFIT attempts to solve problems of the form:min sum {(FUN(X,XDATA)-YDATA).^2} where X, XDATA, YDATA and the valuesX returned by FUN can be vectors ormatrices.X=LSQCURVEFIT(FUN,X0,XDATA,YDATA) starts at X0 and finds coefficients Xto best fit the nonlinear functions in FUN to the data YDATA (in theleast-squares sense). FUN accepts inputs X and XDATA and returns avector (or matrix) of function values F, where F is the same size asYDATA, evaluated at X and XDATA. NOTE: FUN should returnFUN(X,XDATA)and not the sum-of-squares sum((FUN(X,XDATA)-YDATA).^2).((FUN(X,XDATA)-YDATA) is squared and summed implicitly in thealgorithm.)X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB) defines a set of lower andupper bounds on the design variables, X, so that the solution is in therange LB <= X <= UB. Use empty matrices for LB and UB if no boundsexist. Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf ifX(i) is unbounded above.X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB,OPTIONS) minimizes with thedefault parameters replaced by values in the structure OPTIONS, anargument created with the OPTIMSET function. See OPTIMSET for details.Used options are Display, TolX, TolFun, DerivativeCheck, Diagnostics,FunValCheck, Jacobian, JacobMult, JacobPattern, LineSearchType,LevenbergMarquardt, MaxFunEvals, MaxIter, DiffMinChange andDiffMaxChange, LargeScale, MaxPCGIter, PrecondBandWidth, TolPCG,OutputFcn, and TypicalX. Use the Jacobian option to specify that FUNalso returns a second output argument J that is the Jacobian matrix atthe point X. If FUN returns a vector F of m components when X has length n, then J is an m-by-n matrix where J(i,j) is the partialderivative of F(i) with respect to x(j). (Note that the Jacobian J isthe transpose of the gradient of F.)[X,RESNORM]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) returns the valueof thesquared 2-norm of the residual at X: sum {(FUN(X,XDATA)-YDATA).^2}.[X,RESNORM,RESIDUAL]=LSQCURVEFIT(FUN,X0,...) returns the value of residual,FUN(X,XDATA)-YDATA, at the solution X.[X,RESNORM,RESIDUAL,EXITFLAG]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) returnsan EXITFLAG that describes the exit condition of LSQCURVEFIT. Possiblevalues of EXITFLAG and the corresponding exit conditions are1 LSQCURVEFIT converged to a solution X.2 Change in X smaller than the specified tolerance.3 Change in the residual smaller than the specified tolerance.4 Magnitude of search direction smaller than the specified tolerance.0 Maximum number of function evaluations or of iterations reached.-1 Algorithm terminated by the output function.-2 Bounds are inconsistent.-4 Line search cannot sufficiently decrease the residual alongthecurrent search direction.[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=LSQCURVEFIT(FUN,X0,XDATA,YDATA ,...)returns a structure OUTPUT with the number of iterations taken inOUTPUT.iterations, the number of function evaluations inOUTPUT.funcCount,the algorithm used in OUTPUT.algorithm, the number of CG iterations (ifused) in OUTPUT.cgiterations, the first-order optimality (if used)inOUTPUT.firstorderopt, and the exit message in OUTPUT.message.[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA]=LSQCURVEFIT(FUN,X0,XDAT A,YDATA,...)returns the set of Lagrangian multipliers, LAMBDA, at the solution:LAMBDA.lower for LB and LAMBDA.upper for UB.[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN]=LSQCURVEFIT(FU N,X0,XDATA,YDATA,...)returns the Jacobian of FUN at X.ExamplesFUN can be specified using @:xdata = [5;4;6]; % example xdataydata = 3*sin([5;4;6])+6; % example ydatax = lsqcurvefit(@myfun, [2 7], xdata, ydata)where myfun is a MATLAB function such as:function F = myfun(x,xdata)F = x(1)*sin(xdata)+x(2);FUN can also be an anonymous function:x = lsqcurvefit(@(x,xdata) x(1)*sin(xdata)+x(2),[2 7],xdata,ydata)If FUN is parameterized, you can use anonymous functions to capture theproblem-dependent parameters. Suppose you want to solve the curve-fittingproblem given in the function myfun, which is parameterized by its secondargument c. Here myfun is an M-file function such asfunction F = myfun(x,xdata,c)F = x(1)*exp(c*xdata)+x(2);To solve the curve-fitting problem for a specific value of c, first assignthe value to c. Then create a two-argument anonymous function that capturesthat value of c and calls myfun with three arguments. Finally, pass thisanonymous function to LSQCURVEFIT:xdata = [3; 1; 4]; % example xdataydata = 6*exp(-1.5*xdata)+3; % example ydatac = -1.5; % define parameterx = lsqcurvefit(@(x,xdata) myfun(x,xdata,c),[5;1],xdata,ydata) See also optimset, lsqnonlin, fsolve, @, inline.Reference page in Help browserdoc lsqcurvefit>>第二题:1 三次线性拟合clear allx=0:0.5:5;y=x.^3-6*x.^2+5*x-3;y1=y;for i=1:length(y)y1(i)=y1(i)+rand;enda=polyfit(x,y1,3);b=polyval(a,x);plot(x,y,'*',x,b),aa =1.0121 -6.1033 5.1933 -2.4782② 二次线性拟合clear allx=0:0.5:20;y=x.^3-6*x.^2+5*x-3;y1=y;for i=1:length(y)y1(i)=y1(i)+rand;enda=polyfit(x,y1,2);b=polyval(a,x);plot(x,y,'*',x,b),aa =23.9982 -232.0179 367.9756③ 四次线性拟合clear allx=0:0.5:20;y=x.^3-6*x.^2+5*x-3;y1=y;for j=1:length(y)y1(j)=y1(j)+rand;enda=polyfit(x,y1,4);b=polyval(a,x);plot(x,y,'*',x,b),aa =-0.0001 1.0038 -6.0561 5.2890 -2.8249 >>第三题:1 拟合曲线为:f(x)=定义函数:function f=fun(a,x)f=a(1)-(a(1)-a(2))*exp(-a(3)*x);主程序:clear allclcx=[0.5 1 2 3 4 5 7 9];y=[6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63];a0=[1 1 1];a=lsqcurvefit('fun',a0,x,y);y1=a(1)-(a(1)-a(2))*exp(-a(3)*x);plot(x,y,'r*',x,y1,'b')V1=a(2)tei=1/a(3)Optimization terminated: relative function value changing by less than OPTIONS.TolFun.。
matlab最小二乘法曲线拟合
matlab最小二乘法曲线拟合最小二乘法是一种常用的曲线拟合方法,它通过最小化实际观测值与拟合曲线之间的平方误差来确定最佳拟合曲线的参数。
给定一组实际观测数据点(xi, yi),我们希望找到一个拟合曲线y=f(x;θ),其中θ表示曲线的参数。
最小二乘法的目标是使误差的平方和最小化,即使得下述损失函数最小化:L(θ) = ∑(yi - f(xi;θ))^2其中,∑表示求和运算,L(θ)是损失函数,yi是第i个观测数据点的输出值,f(xi;θ)是根据参数θ计算得到的拟合曲线在第i个观测点的预测值。
为了找到最佳的参数θ,我们通过最小化损失函数来求解优化问题:minimize L(θ)这个问题可以通过求解等式∂L/∂θ = 0 来得到最优参数θ的闭式解。
具体的求解方法,可以通过矩阵和向量的运算来实现。
在Matlab中,可以使用“polyfit”函数进行最小二乘法的曲线拟合。
该函数可以拟合出一条多项式曲线,通过指定最佳拟合的次数,即多项式的阶数。
拟合结果包括最佳参数和拟合误差等信息。
使用方法如下:```% 输入观测数据x = [x1, x2, x3, ...]';y = [y1, y2, y3, ...]';% 拟合曲线order = 1; % 最佳拟合的次数(如一次线性拟合)p = polyfit(x, y, order);% 最佳参数coefficients = p;% 拟合曲线curve = polyval(p, x);% 绘制拟合曲线和观测数据plot(x, y, 'o', x, curve)```这样,就可以使用Matlab的最小二乘法曲线拟合方法来得到最佳的拟合曲线。
Matlab最小二乘法曲线拟合
最小二乘法在曲线拟合中比较普遍。
拟合的模型主要有1.直线型2.多项式型3.分数函数型4.指数函数型5.对数线性型6.高斯函数型......一般对于LS问题,通常利用反斜杠运算“\”、fminsearch或优化工具箱提供的极小化函数求解。
在Matlab中,曲线拟合工具箱也提供了曲线拟合的图形界面操作。
在命令提示符后键入:cftool,即可根据数据,选择适当的拟合模型。
“\”命令1.假设要拟合的多项式是:y=a+b*x+c*x^2.首先建立设计矩阵X:X=[ones(size(x)) x x^2];执行:para=X\ypara中包含了三个参数:para(1)=a;para(2)=b;para(3)=c;这种方法对于系数是线性的模型也适应。
2.假设要拟合:y=a+b*exp(x)+cx*exp(x^2)设计矩阵X为X=[ones(size(x)) exp(x) x.*exp(x.^2)];para=X\y3.多重回归(乘积回归)设要拟合:y=a+b*x+c*t,其中x和t是预测变量,y是响应变量。
设计矩阵为X=[ones(size(x)) x t] %注意x,t大小相等!para=X\ypolyfit函数polyfit函数不需要输入设计矩阵,在参数估计中,polyfit会根据输入的数据生成设计矩阵。
1.假设要拟合的多项式是:y=a+b*x+c*x^2p=polyfit(x,y,2)然后可以使用polyval在t处预测:y_hat=polyval(p,t)polyfit函数可以给出置信区间。
[p S]=polyfit(x,y,2) %S中包含了标准差[y_fit,delta] = polyval(p,t,S) %按照拟合模型在t处预测在每个t处的95%CI为:(y_fit-1.96*delta, y_fit+1.96*delta)2.指数模型也适应假设要拟合:y = a+b*exp(x)+c*exp(x.?2)p=polyfit(x,log(y),2)fminsearch函数fminsearch是优化工具箱的极小化函数。
最小二乘法曲线拟合的Matlab程序
方便大家使用的最小二乘法曲线拟合的Matlab程序非常方便用户使用,直接按提示操作即可;这里我演示一个例子:(红色部分为用户输入部分,其余为程序运行的结果,结果图为Untitled.fig,Untitled2.fig) 请以向量的形式输入x,y.x=[1,2,3,4]y=[3,4,5,6]通过下面的交互式图形,你可以事先估计一下你要拟合的多项式的阶数,方便下面的计算.polytool()是交互式函数,在图形上方[Degree]框中输入阶数,右击左下角的[Export]输出图形回车打开polytool交互式界面回车继续进行拟合输入多项式拟合的阶数m = 4Warning: Polynomial is not unique; degree >= number of data points. > In polyfit at 72In zxecf at 64输出多项式的各项系数a = 0.0200000000000001a = -0.2000000000000008a = 0.7000000000000022a = 0.0000000000000000a = 2.4799999999999973输出多项式的有关信息 SR: [4x5 double]df: 0normr: 2.3915e-015Warning: Zero degrees of freedom implies infinite error bounds.> In polyval at 104In polyconf at 92In zxecf at 69观测数据拟合数据x y yh1.0000 3.0000 3.00002.0000 4.0000 4.00003 5 54.0000 6.0000 6.0000剩余平方和 Q = 0.000000标准误差 Sigma = 0.000000相关指数 RR = 1.000000请输入你所需要拟合的数据点,若没有请按回车键结束程序.输入插值点x0 = 3输出插值点拟合函数值 y0 = 5.0000>>结果:untitled.figuntitled2.fig一些matlab优化算法代码的分享代码的目录如下:欢迎讨论1.约束优化问题:minRosen(Rosen梯度法求解约束多维函数的极值)(算法还有bug) minPF(外点罚函数法解线性等式约束)minGeneralPF(外点罚函数法解一般等式约束)minNF(内点罚函数法)minMixFun(混合罚函数法)minJSMixFun(混合罚函数加速法)minFactor(乘子法)minconPS(坐标轮换法)(算法还有bug)minconSimpSearch(复合形法)2.非线性最小二乘优化问题minMGN(修正G-N法)3.线性规划:CmpSimpleMthd(完整单纯形法)4.整数规划(含0-1规划)DividePlane(割平面法)ZeroOneprog(枚举法)5.二次规划QuadLagR(拉格朗日法)ActivedeSet(起作用集法)6.辅助函数(在一些函数中会调用)minNT(牛顿法求多元函数的极值)Funval(求目标函数的值)minMNT(修正的牛顿法求多元函数极值)minHJ(黄金分割法求一维函数的极值)7.高级优化算法1)粒子群优化算法(求解无约束优化问题)1>PSO(基本粒子群算法)2>YSPSO(待压缩因子的粒子群算法)3>LinWPSO(线性递减权重粒子群优化算法)4>SAPSO(自适应权重粒子群优化算法)5>RandWSPO(随机权重粒子群优化算法)6>LnCPSO(同步变化的学习因子)7>AsyLnCPSO(异步变化的学习因子)(算法还有bug)8>SecPSO(用二阶粒子群优化算法求解无约束优化问题)9>SecVibratPSO(用二阶振荡粒子群优化算法求解五约束优化问题)10>CLSPSO(用混沌群粒子优化算法求解无约束优化问题)11>SelPSO(基于选择的粒子群优化算法)12>BreedPSO(基于交叉遗传的粒子群优化算法)13>SimuAPSO(基于模拟退火的粒子群优化算法)2)遗传算法1>myGA(基本遗传算法解决一维约束规划问题)2>SBOGA(顺序选择遗传算法求解一维无约束优化问题)3>NormFitGA(动态线性标定适应值的遗传算法求解一维无约束优化问题)4>GMGA(大变异遗传算法求解一维无约束优化问题)5>AdapGA(自适应遗传算法求解一维无约束优化问题)6>DblGEGA(双切点遗传算法求解一维无约束优化问题)7>MMAdapGA(多变异位自适应遗传算法求解一维无约束优化问题)自己编写的马尔科夫链程序A 代表一组数据序列一维数组本程序的操作对象也是如此t=length(A); % 计算序列“A”的总状态数B=unique(A); % 序列“A”的独立状态数顺序,“E”E=sort(B,'ascend');a=0;b=0;c=0;d=0;for j=1:1:ttLocalization=find(A==E(j)); % 序列“A”中找到其独立状态“E”的位置for i=1:1:length(Localization)if Localization(i)+1>tbreak; % 范围限定elseif A(Localization(i)+1)== E(1)a=a+1;elseif A(Localization(i)+1)== E(2)b=b+1;elseif A(Localization(i)+1)== E(3)c=c+1;% 依此类推,取决于独立状态“E”的个数elsed=d+1;endendT(j,1:tt)=[a,b,c,d]; % “T”为占位矩阵endTT=T;for u=2:1:ttTT(u,:)= T(u,:)- T(u-1,:);endTT; % 至此,得到转移频数矩阵Y=sum(TT,2);for uu=1:1:ttTR(uu,:)= TT(uu,:)./Y(uu,1);endTR % 最终得到马尔科夫转移频率/概率矩阵% 观测序列马尔科夫性质的检验:N=numel(TT);uuu=1;Col=sum(TT,2); % 对列求和Row=sum(TT,1); % 对行求和Total=sum(Row); % 频数总和for i=1:1:ttfor j=1:1:ttxx(uuu,1)=sum((TT(i,j)-(Row(i)*Col(j))./Total).^2./( (Row(i)*Col(j)). /Total));uuu=uuu+1; % 计算统计量x2endendxx=sum(xx)。
最小二乘法拟合matlab
最小二乘法拟合matlab
最小二乘法拟合MATLAB
最小二乘法是一种有效地估计未知参数值的统计学方法,它假定误差服从正态分布,然后进行极大似然估计。
下面我们就来介绍一下如何使用MATLAB来拟合最小二乘法。
1.第一步:绘制出要拟合的数据,这里我们绘制出了一个简单的抛物线数据:
x=[-3 -2 -1 0 1 2 3];
y=[6 3 1 0 -2 -4 -7];
plot(x,y);
2.第二步:根据你要拟合的函数,构建出你所要拟合的模型。
这里,我们想拟合一条抛物线:y=ax2+bx+c ;
3.第三步:定义拟合函数:
fun=@(x,xdata)x(1)*xdata.^2+x(2)*xdata+x(3);
4.第四步:调用最小二乘法函数:
[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcur vefit(fun,[1 1 1],x,y);
现在你已经可以看到拟合函数的参数了:
x的值为[1.7, 0.3, -1.5],
而拟合函数为: y=1.7x2+0.3x-1.5
因此,使用MATLAB调用最小二乘法可以很方便地拟合出任意复
杂的函数,并且可以得到准确的参数值。
Matlab最小二乘法曲线拟合的应用实例
MATLAB机械工程最小二乘法曲线拟合的应用实例班级:姓名:学号:指导教师:一,实验目的通过Matlab上机编程,掌握利用Matlab软件进行数据拟合分析及数据可视化方法二,实验内容1.有一组风机叶片的耐磨实验数据,如下表所示,其中X为使用时间,单位为小时h,Y为磨失质量,单位为克g。
要求:对该数据进行合理的最小二乘法数据拟合得下列数据。
x=[10000 11000 12000 13000 14000 15000 16000 170 00 18000 19000 20000 21000 22000 23000];y=[24.0 26.5 29.8 32.4 34.7 37.7 41.1 42.8 44.6 47.3 6 5.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,6 5.8,87.5,137.8,174.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;endQ=1-chi2cdf(chi2,freedom); %判断拟合良好度clf,shgsubplot(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 onerrorbar(x,YE{nod},D{nod},‘r’);hold offtitle(‘较适当阶次的拟合’)text(10000,150.0,[‘chi2=’num2str(chi2(nod))‘~’int2str(freedom(nod))])text(10000,140.0,[‘freedom=’int2str(freedom(nod))]) text(20000,40.0,[‘Q=’num2str(Q(nod))‘~0.5’])disp(‘’)disp(‘拟合多项式系数’),disp(A{nod})disp(‘拟合系数的离差’),disp(DA{nod})运行结果分为两个阶段,第一阶段先判断拟合度,第二阶段根据拟合度,选择合适的拟合阶次,再绘出拟合结果。
matlab 多点利用最小二乘法拟合二次函数方程
Matlab是一种用于数学计算和工程䇹算的高级语言和交互式环境。
在Matlab中,利用最小二乘法来拟合二次函数方程是一种常见的数据分析方法,可以通过拟合得到二次函数的系数,从而更好地理解和分析实际问题中的数据。
1. 理论基础最小二乘法是一种数学优化方法,用于寻找真实数据与拟合函数之间的最小误差。
在拟合二次函数方程时,我们可以将拟合方程写成如下形式:y = a*x^2 + b*x + c其中,a、b、c分别为二次函数的系数,x和y分别为自变量和因变量。
2. Matlab中的多点利用最小二乘法在Matlab中,可以使用polyfit函数来实现对多点数据进行二次函数拟合。
其基本语法为:p = polyfit(x, y, n)其中,x和y分别为输入的数据点,n为二次函数的次数。
3. 示例代码下面给出一个简单的示例代码来演示如何在Matlab中进行多点利用最小二乘法拟合二次函数方程:```Matlabx = [1, 2, 3, 4, 5];y = [3.1, 4.9, 7.2, 9.8, 12.5];p = polyfit(x, y, 2);```在这个例子中,我们输入了5个数据点,然后利用polyfit函数对这些数据点进行二次函数拟合,得到了二次函数的系数p。
4. 结果分析经过拟合得到的二次函数系数p为:p = [0.1, 0.2, 3]这意味着拟合得到的二次函数方程为:y = 0.1*x^2 + 0.2*x + 3通过这个拟合方程,我们可以更好地理解和分析实际数据的趋势和规律。
5. 需要注意的问题在利用最小二乘法拟合二次函数方程时,需要注意以下几个问题:- 数据的选择:数据点的选择对拟合结果会有很大的影响,需要根据实际问题合理选择数据点。
- 拟合精度:拟合得到的二次函数方程的精度取决于数据的质量和数量,需要谨慎选择拟合的次数。
利用最小二乘法在Matlab中拟合二次函数方程是一种常见且有效的数据分析方法,通过对实际数据进行拟合,可以更好地理解和分析数据规律。
基于matlab的最小二乘法应用
fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy), E2=sqrt(sum(fy2))
X=0:1:55;
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('数据点(xi,yi)和拟合曲线y=f(x)的图形')
方案一:设想 具有指数形式
(4-1)
在求解参数a和b时,为了避免求解一个非线性方程组,对上式两边取对数方程变为:
(4-2)
引入新的变量 ,并记 上式变为 ,此时的问题就转化为求形如 的最小二乘解。运用matlab语言编写计算和画图程序,程序一见附录部分[2]。运算的结果:a=5.2151,b=-7.4962。最大偏差 ,均方误差 。故拟合的曲线为:
(2-7)
其中
(2-8)
它的均方误差为: (2-9)
最大偏差为: (2-10)
3
在化学反应中,由实验测得分解物浓度与时间的关系如下表2所示
表2 浓度(y)与时间( x )的关系实验数据表
x
0
5
10
15
20
25
y
0
1.27
2.16
2.86
3.44
3.87
x
30
35
40
45
50
曲线拟合的最小二乘法matlab举例
曲线拟合的最小二乘法学院:光电信息学院 姓名:赵海峰 学号:200820501001一、曲线拟合的最小二乘法原理:由已知的离散数据点选择与实验点误差最小的曲线)(...)()()(1100x a x a x a x S n n ϕϕϕ+++=称为曲线拟合的最小二乘法。
若记),()()(),(0i k i j mi i k j x x x ϕϕωϕϕ∑==k i k i mi i k d x x f x f ≡=∑=)()()(),(0ϕωϕ上式可改写为),...,1,0(;),(n k d a k j noj j k -=∑=ϕϕ这个方程成为法方程,可写成距阵形式d Ga =其中,),...,,(,),...,,(1010T n T n d d d d a a a a ==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=),(),(),()(),(),(),(),(),(101110101000n n n n n n G ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ 。
它的平方误差为:.)]()([)(||||2022i i mi i x f x S x -=∑=ωδ二、数值实例:下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度数据表,按照数据找出任意次曲线拟合方程和它的图像。
下面应用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]; a1=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).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%四、数值结果:不同次数多项式拟和误差平方和为:r1 = 67.6659r2 = 20.1060r3 = 3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
matlab 最小二乘法拟合椭圆
matlab 最小二乘法拟合椭圆在MATLAB中,可以使用最小二乘法对一组数据进行椭圆拟合。
最小二乘法是一种常见的数值拟合方法,通过最小化实际数据点与拟合曲线之间的差异来确定最佳拟合参数。
首先,将椭圆的方程表示为:(x - h)^2 / a^2 + (y - k)^2 / b^2 = 1其中(h, k)是椭圆的中心坐标,a和b是椭圆的半长轴和半短轴长度。
令数据点的坐标为(xi, yi),通过最小化以下误差函数来拟合椭圆:F = sum(((xi - h)^2 / a^2 + (yi - k)^2 / b^2) - 1)^2其中,求和遍历所有数据点。
为了找到最佳的拟合参数h、k、a和b,可以使用MATLAB中的最小二乘法拟合函数如lsqcurvefit。
以下是使用最小二乘法进行椭圆拟合的MATLAB代码示例:```Matlab% 假设有一组包含椭圆上的数据点的二维矩阵data,每行包含一个点的坐标(xi, yi)% 定义误差函数fun = @(params, x) ((x(:, 1) - params(1)).^2 ./ params(3)^2 + (x(:, 2) - params(2)).^2 ./ params(4)^2 - 1).^2;% 初始化参数的初始猜测值params0 = [0, 0, 1, 1];% 使用最小二乘法进行拟合params = lsqcurvefit(fun, params0, data(:, 1), data(:, 2));% 提取拟合的椭圆参数h = params(1); % 中心坐标xk = params(2); % 中心坐标ya = params(3); % 半长轴长度b = params(4); % 半短轴长度% 绘制原始数据点和拟合的椭圆figure;plot(data(:, 1), data(:, 2), 'bo'); % 原始数据点hold on;theta = linspace(0, 2*pi, 100);x = h + a*cos(theta); % x坐标y = k + b*sin(theta); % y坐标plot(x, y, 'r-', 'LineWidth', 2); % 拟合的椭圆axis equal;xlabel('x');ylabel('y');title('椭圆拟合');legend('数据点', '拟合椭圆');```在以上代码中,首先定义了误差函数fun,该函数计算数据点与拟合椭圆之间的差异。
曲线拟合的最小二乘法的Matlab程序实现
《曲线拟合的最小二乘法》目录1 实验要求 (1)1.1 问题提出 (1)1.2 要求 (1)1.3 目的和意义 (1)2 问题求解 (2)2.1 基本原理 (2)2.1.1 最小二乘法简介 (2)2.1.2 求解思路 (2)2.1.3 常用的多项式拟合 (4)2.2 求解步骤 (5)2.2.1 绘制散点图 (5)2.1.2 进行三次拟合 (6)2.1.3 三次拟合的误差分析 (7)2.1.4 进行五次拟合 (8)2.1.5 五次拟合的误差分析 (10)2.1.6 拟合效果比较 (11)3 参考文献 (13)1 实验要求1.1 问题提出从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量 y 与时间t 的拟合曲线。
1.2 要求(1)用最小二乘法进行曲线拟合;(2)近似解析表达式为φ(t)=a1t+a2t2+a3t3;(3)打印出拟合函数φ(t),并打印出φ(t j) 与y(t j) 的误差,j=1,2,⋯,12 ;(4)另外选取一个近似表达式,尝试拟合效果的比较;(5) * 绘制出曲线拟合图﹡。
1.3 目的和意义(1)掌握曲线拟合的最小二乘法;(2)最小二乘法亦可用于解超定线代数方程组;(3)探索拟合函数的选择与拟合精度间的关系。
2 问题求解2.1 基本原理2.1.1 最小二乘法简介对原始数据做适当变换或其他处理,发现其中隐藏的数学规律往往是数学建模中非常重要的一步。
通过对实验数据的分析,求出自变量、因变量之间近似的函数关系,进而对其进行进一步分析,而最小二乘法是对数据进行拟合常用的方法。
最小二乘法是研究观测数据的主要方法之一,若已知两变量满足线性关系y =ax +b ,对这两个变量进行了 n 次观测,从而得到了 n 对观测数据(x 1, y 1),(x 2, y 2)......(x n , y n ),最小二乘法提供了一个方法,基本思想是找到离这 n 组观测数据最近的直线。
用MatLab画图(最小二乘法做曲线拟合)
---------------------------------------------------------------最新资料推荐------------------------------------------------------ 用MatLab画图(最小二乘法做曲线拟合) 用 MatLab 画图(最小二乘法做曲线拟合) 帮朋友利用实验数据画图时,发现 MatLab 的确是画图的好工具,用它画的图比Excel光滑、精确。
利用一组数据要计算出这组数据对应的函数表达式从而得到相应图像,MatLab 的程序如下:x=[1 5 10 20 30 40 60 80] y=[15. 4 33. 9 42. 2 50. 556 62. 7 72 81. 1] plot(x, y, ‘ r*’ ) ; legend(‘ 实验数据(xi, yi) ‘ ) xlabel(‘ x’ ) , ylabel(‘ y’ ) , title(‘ 数据点(xi, yi) 的散点图’ ) syms a1 a2 a3 x=[15 10 20 30 40 60 80]; fi=a1. *x. +a2. *x+a3 y=[15. 4 33.9 42. 2 50. 5 56 62. 7 72 81. 1] fi =[a1+a2+a3,25*a1+5*a2+a2+(400*a1+20*a2+a3-101/2) +(900*a1+30*a2+a3-56) +(1600*a1+40*a2+a3-627/10) +(3600*a1+60*a2+a3-72)+(6400*a1+80*a2+a3-811/10) ; Ja1=diff(J, a1) ;Ja2=diff(J, a2) ; Ja3=diff(J, a3) ; Ja11=simple(Ja1) ,Ja21=simple(Ja2) , Ja31=simple(Ja3) A=[114921252, 1656252, 26052; 1656252, 26052, 492; 26052, 492, 16]; B=[9542429/5, 166129/5, 4138/5]; C=B/A, f=poly2sym(C) xi=[1 5 10 20 3040 60 80] ; y=[15. 4 33. 9 42. 2 50. 5 56 62. 7 72 81. 1]; n=length(xi) ; f=-0. 0086. *xi. +1. 3876. *xi+23. 1078;1 / 6x=1: 1/10: 80; F=-0. 0086. *x. +1. 3876. *x+23. 1078; fy=abs(f-y) ; fy2=fy. ; 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) ‘ , ‘ 拟合曲线f(x) = -0. 0086x +1. 3876x+23. 1078’ ) , xlabel(‘ x’ ) , ylabel(‘ y’ ) , title(‘ 实验数据点(xi, yi) 及拟合曲线f(x) ‘ ) 下图是程序运行后得到的:Su7Tw8VxaW ybXAcZB d#Cf!Eg%FhGj*Ik(Jl-Kn+M o0Np2Or3Ps4R t6Sv7Tw8VxaWzbXAcZBe#Cf! Eg%Fi Gj*Ik)Jl-Kn+Mo1Np2Or3Qs4Rt6Sv7 Uw8Vx aWzbYAc ZBe#Df!Eg%FiHj*Ik) Jm-Kn +Mo1Nq2Or3Qs 5Rt6Sv7Uw9VxaWzbYAdZBe#D f$Eg%F iHj(I k) Jm-Ln+Mo1Nq2Pr3Qs5Ru6S v7Uw9V yaWzbY AdZCe#Df$Eh%FiHj(Il) Jm-Ln0Mo1Nq2Pr4 Qs5Ru6Tv8Uw9VyaXzbYAdZCe !Df$Eh %GiHj (Il) Km-Ln0Mp1Nq2Pr4Qt5Ru 6Tv8U x9VyaXz cYAdZCe! Dg$Eh%Gi*Hj(Il) Km+Ln0M p1Oq2P r4Qt5Su6Tv8Ux9WyaXzcYBdZ Ce!Dg$Fh%Gi* Hk(Il) Km+Lo0Mp1Oq3Pr4Qt5 Su7Tv8Ux9Wyb XzcYBd#Ce!Dg$FhGi*Hk(Jl ) Km+L o0Np1Oq 3Ps4Rt 5Su7Tw8Ux9WybXAcY Bd#Cf!Dg$FhGj*Hk(Jl-Km+Lo0Np2Oq3Ps4 Rt6Su7Tw8Vx9 WybXAcZBd#Cf!Eg$FhGj*Ik (Jl-Kn +Lo0Np2Or3Ps4Rt6Sv7Tw8VxaWybXA cZBe#Cf!Eg%F hGj*Ik) Jl-K n+Mo0Np2Or3Q s4Rt6Sv 7Uw8V xaWzbXAcZBe# D f! Eg%FiGj* Ik) Jm- Kn+Mo1 Nq2Or3Qs5Rt6Sv7Uw9VxaWzb YAcZBe#Df$Eg %FiHj*Ik)---------------------------------------------------------------最新资料推荐------------------------------------------------------ Jm-Ln+Mo1Nq2Pr 3Qs5Ru 6Sv7Uw 9VyaWzbYAdZBe#Df$Eh%FiH j(Ik) J m-Ln0M o1Nq2Pr4Qs5Ru6Tv7Uw9VyaX zbYAdZC e#Df$ Eh%GiHj(Il) Jm-Ln0Mp1Nq2 Pr4Qt5Ru6Tv8 Uw9VyaXzcYAdZCe!Df$Eh%Gi *Hj(Il) Km+Ln 0Mp1Oq2Pr4Qt5Su6Tv8Ux9Vy aXzcYB dZCe!D g$Eh%Gi*Hk(Il) Km+Lo0Mp1O q3Pr4Qt5Su7Tv8Ux9WyaXzc Y Bd#Ce!Dg$Fh %Gi*Hk( Jl) Km +Lo0Np1Oq3Ps 4 Qt5Su7Tw8Ux 9WybXzcYBd#C f!Dg$FhGi*H k (Jl-Km+Lo0N p2Oq3Ps4Rt5S u7Tw8Vx9WybX AcYBd#Cf! Eg$ FhGj*Ik (Jl- Kn+Lo0Np2Or3 P s4Rt6Su7Tw8 VxaWybXA cZBd #Cf!Eg%FhGj * Ik) Jl-Kn+Mo 0Np2Or3Qs4Rt 6Sv7Tw8VxaWz bXAcZBe#Cf!E g%FiGj*Ik) J m-Kn+Mo1Np2O r 3Qs5Rt6Sv7U w8VxaWzbYAcZ Be#Df! Eg%Fi H j*Ik) Jm-Ln+ Mo1Nq2O r3Qs5 Ru6Sv7Uw9Vxa W zbYAdZBe#Df $Eh%Fi Hj(Ik ) Jm-Ln0Mo1Nq 2Pr3Qs5Ru6Tv 7Uw9Vya WzbYA dZCe#Df$Eh%G iHj(Il) Jm-L n0Mp1Nq2Pr4Q s5Ru6Tv8Uw9V ya XzbYAdZCe! Df$Eh%Gi*Hj( Il) Km-Ln0Mp1 Oq 2Pr4 Qt5Ru6Tv8Ux9Vy aXz cYAdZCe!Dg$E h%G i*Hk(Il) K m+Ln0Mp1O q3P r4Qt5Su6Tv8U x9WyaXzcYBd# Ce!Dg$Fh% Gi* Hk(Jl) Km+Lo0 Mp1Oq3Ps4Qt5 Su7Tv8Ux9Wyb XzcYBd#Cf! Dg $F hGi*Hk(Jl -Km+Lo0N p1Oq 3Ps4Rt5Su7Tw 8U x9WybXAcYB d#Cf!Eg$F hG j*Hk(Jl-Kn+Lo0Np2O q3Ps4Rt 6Su7Tw8Vx9WybXAcZBd#Cf!E g%FhGj*Ik(J l-Kn+Mo0Np2O r3Ps4Rt6Sv7T w8Vxa WzbXAcZ Be#Cf! Eg%Fi Gj*Ik) Jl-Kn+ Mo1Np2Or3Qs4 Rt6Sv7Uw8VxaWzbYAcZBe#Df !Eg%FiHj*Ik )3 / 6Jm-Kn+Mo1Nq2Or3Qs5Rt6Sv 7Uw9Vx aWzbYA dZBe#Df$Eg%FiHj(Ik) Jm-L n+Mo1Nq2Pr3Q s5Ru6Sv7Uw9VyaWzbYAdZCe# Df$Eh %FiHj( Il) Jm-Ln0Mo1Nq2Pr4Qs5Ru6 Tv8Uw9VyaXzb YAdZCe!Df$Eh%GiHj(Il) Km -Ln0Mp 1Nq2Pr 4Qt5Ru6Tv8Ux9VyaXzcYAdZC e!Dg$E h%Gi*H j(Il) Km+Ln0Mp1Oq2Pr4Qt5S u6Tv8U x9WyaX zcYBdZCe! Dg$Fh%Gi*Hk(Il) Km+Lo0Mp1Oq3 Pr4Qt5Su7Tv8Ux9WybXzcYBd #Ce!D g$FhGi *Hk(Jl ) Km+Lo0Np1Oq3Ps4R t5Su7Tw8Ux9T v7Uw9VyaXzbYAdZCe#Df$Eh% GiHj( Il) Jm- Ln0Mp1Nq2Pr4Qs5Ru6Tv8Uw9 VyaXzcY AdZCe !Df$Eh%Gi*Hj(Il) Km-Ln0Mp 1Oq2Pr 4Qt5Ru 6Tv8Ux9VyaXzcYBdZCe! Dg$E h%Gi*Hk(Il) K m+Ln0Mp1Oq3Pr4Qt5Su6Tv8U x9WyaX zcYBd# Ce!Dg$Fh%Gi*Hk(Jl) Km+Lo0 Mp1Oq3P s4Qt5 Su7Tw8Ux9Wyb X zcYBd#Cf!Dg $FhGi*Hk(Jl -Km+Lo0Np1Oq3Ps4Rt5Su7Tw 8Vx9Wy bXAcYB d#Cf! Eg$FhGj*Hk(Jl-Kn+L o0Np2O q3Ps4R t6Su7Tw8VxaWybXAcZBd#Cf! Eg%Fh Gj*Ik( Jl-Kn+Mo0Np2Or3Ps4Rt6Sv7 Tw8Vxa WzbXAc ZBe#Cf!Eg%FiGj*Ik) Jl-Kn +Mo1Np2Or3Qs 5Rt6Sv7Uw8Vx a WzbYAcZBe#D f! Eg%FiHj*Ik) Jm-Kn+Mo1 Nq 2Or3Qs5Ru6 Sv7Uw9Vx aWzb YAdZBe#Df$Eg %F iHj(Ik) Jm -Ln+Mo1N q2Pr 3Qs5Ru6Tv7Uw 9VyaWzbYAdZC e#Df$Eh%FiH j(Il) Jm-Ln0M o1Nq2Pr4Qs5R u6Tv8Uw9VyaX zbYAdZCe! Df$ E h%GiHj(Il) Km-Ln0Mp 1Oq2 Pr4Qt5Ru6Tv8 U x9VyaXzcYAd ZCe!Dg$E h%Gi *Hj(Il) Km+Ln 0Mp1Oq3Pr4Qt 5Su6Tv8Ux9Wy aXzcYBdZCe!D g $Fh%Gi*Hk(I l)---------------------------------------------------------------最新资料推荐------------------------------------------------------ Km+Lo0Mp1O q3Ps4Qt5Su7T v 8Ux9WybXzcY Bd#Ce! D g$Fh Gi*Hk(Jl) Km+ L o0Np1Oq3Ps4 Rt5Su7T w8Ux9 WybXAcYBd#Cf !Dg$FhGj*Hk (Jl-Km+ Lo0Np 2Oq3Ps4Rt6Su 7Tw8Vx9WybXA cZBd#Cf ! Eg$F hGj*Ik(Jl-K n+Mo0Np2Or3P s4Rt6Sv 7Tw8V xaWybXAcZBe# C f!Eg% FhGj*Ik) Jl-K n+Mo 1Np2Or3Qs4Rt 6Sv 7Uw8VxaWz bXAcZBe#D f!E g%FiGj*Ik) J m- Kn+Mo1Nq2O r3Qs5Rt6Sv7U w9VxaWzbYAcZ Be#Df$Eg%Fi Hj*Ik) Jm-Ln+ Mo1Nq2Pr3Qs5 Ru6Sv7Uw9Vya WzbYAdZBe#Df $Eh%FiHj(Il ) Jm -Ln0Mo1Nq 2Pr4Qs5Ru 6Tv 7Uw9VyaXzbYA dZCe#Df$Bd#C f! Eg%FhGj*I k(Jl-Kn+Mo0Np2O r3Qs4Rt 6Sv7Tw8VxaWzbXAcZBe#Cf!E g%FiGj*Ik) J l-Kn+Mo1Np2O r3Qs5Rt6Sv7U w8Vxa WzbYAcZ Be#Df! Eg%Fi Hj*Ik) Jm-Kn+ Mo1Nq2Or3Qs5 Ru6Sv7Uw9VxaWzbYAdZBe#Df $Eg%FiHj(Ik ) Jm-Ln0Mo1Nq2Pr3Qs5Ru6Tv 7Uw9Vy aWzbYA dZCe#Df$Eh%FiHj(Il) Jm-L n0Mp1Nq2Pr4Q s5Ru6Tv8Uw9VyaXzbYAdZCe! Df$Eh %GiHj( Il) Km-Ln0Mp1Oq2Pr4Qt5Ru6 Tv8Ux9VyaXzc YAdZCe!Dg$Eh%Gi*Hj(Il) Km +Ln0Mp 1Oq3Pr 4Qt5Su6Tv8Ux9WyaXzcYBdZC e!Dg$F h%Gi*H k(Jl) Km+Lo0Mp1Oq3Ps4Qt5S u7Ts5R u6Sv7U w9VyaWzbYAdZBe#Df$Eh%Fi Hj(Ik)Jm-Ln0 Mo1Nq2Pr4Qs5Ru6Tv7Uw9Vya XzbYA dZCe#Df $Eh%Gi Hj(Il) Jm-Ln0Mp1N q2Pr4Qt5Ru6T v8Uw9VyaXzcYAdZCe!Df$Eh% Gi*Hj( Il) Km- Ln0Mp1Oq2Pr4Qt5Su6Tv8Ux9 VyaXzcY BdZCe !Dg$Eh%Gi*Hk(Il) Km+Ln0Mp 1Oq3Pr4Qt5Su5 / 67Tv8Ux9WyaXzcYBd#Ce! Dg$F h%Gi*Hk (Jl) K m+Lo0Np1Oq3Ps4Qt5Su7Tw8U x9WybX zcYBd# Cf!Dg$FhGi*Hk(Jl-Km+Lo0 Np2Oq3Ps4Rt5 Su7Tw8Vx9WybXAcYBd#Cf!Eg $FhGj*Hk(Jl -Kn+Lo0Np2Or3Ps4Rt6Su7Tw 8VxaWy bXAcZB d#Cf! Eg%FhG j*Ik(Jl-Kn+M o0Np2O r3Qs4R t6Sr4Qt5Su6Tv8Ux9WyaXzcY BdZCe!Dg$Fh% Gi*Hk(Il) Km+Lo0Mp1Oq3Pr4 Qt5Su7Tv8Ux9 WybXzcYBd#Ce! Dg$FhGi*Hk (Jl) Km +Lo0Np 1Oq3Ps4Rt5Su7Tw8Ux9WybXA cYBd#Cf!Dg$FhGj*Hk(Jl- K m+Lo0Np2Oq3 Ps4Rt6Su7Tw8 Vx9WybXAcZBd # Cf!Eg$FhGj *Ik(Jl- Kn+Lo 0Np2Or3Ps4Rt 6Sv7Tw8VxaWy bXAcZBe#Cf!E g%FhGj*Ik) J l-Kn+Mo0Np2O r3Qs4Rt6Sv7U w8VxaWzbXAcZ B e#Df!Eg%Fi Gj*Ik) J m-Kn+ Mo1Nq2Or3Qs5 R t6Sv7Uw9Vxa WzbYAcZB e#Df $Eg%Ff! Dg$Fh Gi*Hk(Jl-Km +Lo0Np1Oq3Ps 4Rt5Su7Tw8Vx 9W ybXAcYBd#C f!Eg$Fh Gj*H k(Jl-Kn+Lo0N p2Oq3Ps4Rt6S u7Tw8Vxa WybX AcZBd#Cf!Eg% F hGj*Ik(Jl- Kn+Mo0N p2Or3 Ps4Rt6Sv7Tw8 V xaWzbXAcZBe #Cf!Eg% FiGj *Ik) Jl-Kn+Mo 1Np2Or3Qs5Rt 6Sv7Uw8V xaWz bYAcZBe#Df!E g%FiHj*Ik) J m-Kn+Mo1Nq2O r3Qs5Ru6Sv7U w9VxaW zbYAdZBe#Df$Eg %Fi Hj(Ik) Jm-Ln +M o1Nq2Pr3Qs 5Ru6Tv7U w9Vy aWzbYAdZCe#D f$Eh%FiHj(I l) Jm-Ln0Mo1N q2Pr4Qs5Or3P s4Rt6Su7Tw8V xaWybXAcZ Be# Cf! Eg%FhGj* I。
用Matlab进行最小二乘法线性拟合(求传感器非线性误差、灵敏度)
%后面的为注释,红色部分代码需要根据实际情况更改
%最小二乘法线性拟合y=ax+b
x=[,1,,2,,3,,4,,5];%自变量
y=[191,321,442,565,686,819,930,1032,1153,1252];%因变量
xmean=mean(x);ymean=mean(y);
sumx2=(x-xmean)*(x-xmean)';
sumxy=(y-ymean)*(x-xmean)';
a=sumxy/sumx2;%解出直线斜率a(即传感器灵敏度)
b=ymean-a*xmean;%解出直线截距b
z=((a*(x(1,10))+b-(y(1,10)))/(y(1,10)));%“10”是自变量的个数,z为非线性误差(即线性度)
a
b
z
%作图,先把原始数据点用蓝色"十"字描出来
figure
plot(x,y,'+');
hold on
% 用红色绘制拟合出的直线
px=linspace(0,6,50);%(linspace语法(从横坐标负轴起点0画到横坐标正轴终点6,50等分精度)) py=a*px+b;
plot(px,py,'r');
运行结果:
a =
b =
另一种简单一点的方法:
%最小二乘法线性拟合y=ax+b
x=[,1,,2,,3,,4,,5];%自变量
y=[191,321,442,565,686,819,930,1032,1153,1252];%因变量p=polyfit(x,y,1);
p
运行结果:
p =。
最小二乘法曲线拟合的Matlab程序
方便大家使用的最小二乘法曲线拟合的Matlab程序非常方便用户使用,直接按提示操作即可;这里我演示一个例子:(红色部分为用户输入部分,其余为程序运行的结果,结果图为Untitled.fig,Untitled2.fig) 请以向量的形式输入x,y.x=[1,2,3,4]y=[3,4,5,6]通过下面的交互式图形,你可以事先估计一下你要拟合的多项式的阶数,方便下面的计算.polytool()是交互式函数,在图形上方[Degree]框中输入阶数,右击左下角的[Export]输出图形回车打开polytool交互式界面回车继续进行拟合输入多项式拟合的阶数m = 4Warning: Polynomial is not unique; degree >= number of data points. > In polyfit at 72In zxecf at 64输出多项式的各项系数a = 0.0200000000000001a = -0.2000000000000008a = 0.7000000000000022a = 0.0000000000000000a = 2.4799999999999973输出多项式的有关信息 SR: [4x5 double]df: 0normr: 2.3915e-015Warning: Zero degrees of freedom implies infinite error bounds.> In polyval at 104In polyconf at 92In zxecf at 69观测数据拟合数据x y yh1.0000 3.0000 3.00002.0000 4.0000 4.00003 5 54.0000 6.0000 6.0000剩余平方和 Q = 0.000000标准误差 Sigma = 0.000000相关指数 RR = 1.000000请输入你所需要拟合的数据点,若没有请按回车键结束程序.输入插值点x0 = 3输出插值点拟合函数值 y0 = 5.0000>>结果:untitled.figuntitled2.fig一些matlab优化算法代码的分享代码的目录如下:欢迎讨论1.约束优化问题:minRosen(Rosen梯度法求解约束多维函数的极值)(算法还有bug) minPF(外点罚函数法解线性等式约束)minGeneralPF(外点罚函数法解一般等式约束)minNF(内点罚函数法)minMixFun(混合罚函数法)minJSMixFun(混合罚函数加速法)minFactor(乘子法)minconPS(坐标轮换法)(算法还有bug)minconSimpSearch(复合形法)2.非线性最小二乘优化问题minMGN(修正G-N法)3.线性规划:CmpSimpleMthd(完整单纯形法)4.整数规划(含0-1规划)DividePlane(割平面法)ZeroOneprog(枚举法)5.二次规划QuadLagR(拉格朗日法)ActivedeSet(起作用集法)6.辅助函数(在一些函数中会调用)minNT(牛顿法求多元函数的极值)Funval(求目标函数的值)minMNT(修正的牛顿法求多元函数极值)minHJ(黄金分割法求一维函数的极值)7.高级优化算法1)粒子群优化算法(求解无约束优化问题)1>PSO(基本粒子群算法)2>YSPSO(待压缩因子的粒子群算法)3>LinWPSO(线性递减权重粒子群优化算法)4>SAPSO(自适应权重粒子群优化算法)5>RandWSPO(随机权重粒子群优化算法)6>LnCPSO(同步变化的学习因子)7>AsyLnCPSO(异步变化的学习因子)(算法还有bug)8>SecPSO(用二阶粒子群优化算法求解无约束优化问题)9>SecVibratPSO(用二阶振荡粒子群优化算法求解五约束优化问题)10>CLSPSO(用混沌群粒子优化算法求解无约束优化问题)11>SelPSO(基于选择的粒子群优化算法)12>BreedPSO(基于交叉遗传的粒子群优化算法)13>SimuAPSO(基于模拟退火的粒子群优化算法)2)遗传算法1>myGA(基本遗传算法解决一维约束规划问题)2>SBOGA(顺序选择遗传算法求解一维无约束优化问题)3>NormFitGA(动态线性标定适应值的遗传算法求解一维无约束优化问题)4>GMGA(大变异遗传算法求解一维无约束优化问题)5>AdapGA(自适应遗传算法求解一维无约束优化问题)6>DblGEGA(双切点遗传算法求解一维无约束优化问题)7>MMAdapGA(多变异位自适应遗传算法求解一维无约束优化问题)自己编写的马尔科夫链程序A 代表一组数据序列一维数组本程序的操作对象也是如此t=length(A); % 计算序列“A”的总状态数B=unique(A); % 序列“A”的独立状态数顺序,“E”E=sort(B,'ascend');a=0;b=0;c=0;d=0;for j=1:1:ttLocalization=find(A==E(j)); % 序列“A”中找到其独立状态“E”的位置for i=1:1:length(Localization)if Localization(i)+1>tbreak; % 范围限定elseif A(Localization(i)+1)== E(1)a=a+1;elseif A(Localization(i)+1)== E(2)b=b+1;elseif A(Localization(i)+1)== E(3)c=c+1;% 依此类推,取决于独立状态“E”的个数elsed=d+1;endendT(j,1:tt)=[a,b,c,d]; % “T”为占位矩阵endTT=T;for u=2:1:ttTT(u,:)= T(u,:)- T(u-1,:);endTT; % 至此,得到转移频数矩阵Y=sum(TT,2);for uu=1:1:ttTR(uu,:)= TT(uu,:)./Y(uu,1);endTR % 最终得到马尔科夫转移频率/概率矩阵% 观测序列马尔科夫性质的检验:N=numel(TT);uuu=1;Col=sum(TT,2); % 对列求和Row=sum(TT,1); % 对行求和Total=sum(Row); % 频数总和for i=1:1:ttfor j=1:1:ttxx(uuu,1)=sum((TT(i,j)-(Row(i)*Col(j))./Total).^2./( (Row(i)*Col(j)). /Total));uuu=uuu+1; % 计算统计量x2endendxx=sum(xx)。
matlab最小二乘法拟合反应谱代码
一、介绍最小二乘法拟合反应谱最小二乘法是一种常见的数学拟合方法,它可以用来寻找一组数据点的最佳拟合曲线或曲面。
在地震工程中,反应谱是一种用来描述建筑结构在地震作用下的振动情况的重要工具。
使用最小二乘法可以对地震反应谱进行拟合,从而得到结构的振动特性和相应的参数。
二、Matlab中的最小二乘法拟合反应谱代码在Matlab中,有许多函数和工具箱可以用来进行最小二乘法拟合反应谱。
下面是一个示例代码,用来拟合地震反应谱,并得到相关的参数。
```matlab读入反应谱数据T = [0.01 0.02 0.03 0.04 0.05]; 周期Sa = [0.1 0.15 0.2 0.25 0.3]; 加速度反应谱使用最小二乘法拟合反应谱p = polyfit(log(T), log(Sa), 1); 进行对数拟合a = exp(p(2)); 求幂b = p(1); 求斜率输出拟合结果disp(['拟合的函数为:Sa = ' num2str(a) '*T^' num2str(b)]);```在上面的代码中,首先读入了一组反应谱数据,然后使用polyfit函数进行最小二乘法拟合。
接着利用拟合结果得到了拟合的函数,并将其输出。
三、代码解析1. 读入反应谱数据在实际应用中,首先需要将反应谱数据保存在数组中,包括周期T和加速度反应谱Sa。
这些数据可能来自于实测或者地震波形分析得到。
2. 使用最小二乘法拟合反应谱在Matlab中,可以利用polyfit函数对数据进行最小二乘法拟合。
在上面的示例代码中,使用了polyfit(log(T), log(Sa), 1)进行对数拟合。
其中log(T)和log(Sa)是对原始数据取对数,然后利用polyfit进行线性拟合得到拟合参数。
3. 输出拟合结果利用拟合参数得到了拟合的函数,并将其输出。
在本例中,输出结果为Sa = a*T^b,其中a和b分别是拟合得到的系数。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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))])
text(10000,140.0,[‘freedom=’int2str(freedom(nod))])
text(20000,40.0,[‘Q=’num2str(Q(nod))‘~0.5’])
disp(‘’)
disp(‘拟合多项式系数’),disp(A{nod})
disp(‘拟合系数的离差’),disp(DA{nod})
运行结果分为两个阶段,第一阶段先判断拟合度,第二阶段根据拟合度,选择合适的拟合阶次,再绘出拟合结果。
运行结果如下:
>>根据图形选择适当的阶次(请输入数值)4
拟合多项式系数
0.0000 -0.0000 0.0000 -0.1609 504.2295
拟合系数的离差
0.0000 0.0000 0.0000 0.0399 153.3130
四.输出图像。