最小二乘法程序说明及流程图

陆韶琦 3110000441

程序说明：本程序用多项式拟合数据，程序会要求输入需要拟合的次数和数据点的个数，数据文件应该保存在本程序运行时的current folder下，文件取名为“mytext.txt”

程序代码：

%多项式最小二乘法拟合数据

N=input('please put in how many times the power will you overfit:'); M=input('how many couples of statistics are there in the table:');

%读入数据文件

f=fopen('mytxt.txt','r');

S=fscanf(f,'%g',[M 2]);

fclose(f);

S=S';

%显示数据文件，确保正确输入

disp('S(x,y)=');

disp(S);

%建立多项式系数法方程组中间矩阵

C=zeros(N+1,M);

for i=1:N+1

for j=1:M

if S(1,j)==0

C(i,j)=0;

else

C(i,j)=S(1,j).^(i-1);

end

end

end

%建立法方程组

A=C*C';

Y=zeros(M,1);

for i=1:M

Y(i,1)=S(2,i);

end

b=C*Y;

%用列主元高斯消元法接法方程组

A=[A,b];

for i=1:N+1

max=abs(A(i,i));

for j=i+1:N+1

if abs(A(j,i))>max

flag=j;

max=A(j,i);

end

end

for k=i:N+2

B=A(flag,k);

A(flag,k)=A(i,k);

A(i,k)=B;

end

for kh=i+1:N+1

m=-A(kh,i)/A(i,i);

A(kh,i)=0;

for kl=i+1:N+2

A(kh,kl)=A(kh,kl)+m*A(i,kl);

end

end

end

X=zeros(N+1,1);

for i=N+1:-1:1

for j=i-1:-1:1

m=-A(j,i)/A(i,i);

A(j,N+2)=A(j,N+2)+m*A(i,N+2);

end

X(i,1)=A(i,N+2)/A(i,i);

end

disp(X);

%根据系数求得待定曲线

syms x;

expr=0;

for i=1:N+1

expr=expr+X(i,1)*x.^(i-1);

end

%输出得到的曲线表达式

disp(expr);

%计算偏差

bias=zeros(M,1);

for j=1:M

for i=1:N+1

bias(j,1)=bias(j,1)+X(i,1)*S(1,j)^(i-1); end

bias(j,1)=bias(j,1)-S(2,j);

end

%寻找最大偏差

max=abs(bias(1,1));

flag=1;

for i=2:M

if abs(bias(i,1))>max

flag=i;

max=abs(bias(i,1));

end

end

disp('the maximun absoulute value is:'); disp(max);

%计算均方误差

rms=0;

for i=1:M

rms=rms+bias(i,1)^2;

end

rms=sqrt(rms);

disp('the square bias is:');disp(rms);

%制图

a=S(1,1):0.01:S(1,M);

y=subs(expr,x,a);

plot(a,y);

hold on;

grid on;

for i=1:M

x=S(1,i);

y=S(2,i);

plot(x,y,'*');

hold on;

end

运行结果：

表达式中分式难以化

简，但在表达式前给出

了次幂前的四位有效数

字的系数