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

  • 格式:doc
  • 大小:138.00 KB
  • 文档页数:7

下载文档原格式

  / 7
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

用最小二乘法进行多项式拟合(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

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 y i -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552

总共有7个数据点,令m=6

第一步:画出已知数据的的散点图,确定拟合参数n;

x=-1.0:0.5:2.0;y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];

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

end;

再来求矩阵

B=

B=[0 0 0 0];

for j=1:n+1

for i=1:m+1

B(j)=B(j)+y(i)*x(i)^(j-1)

end

end

第三步:写出正规方程,求出a0,,a1…,an. B=B';

a=inv(A)*B;

第四步:画出拟合曲线

x=[-1.0:0.0001:2.0];

z=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3;

plot(x,z)

legend('离散点','y=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3')

title('拟合图')

总程序附下:

x=-1.0:0.5:2.0;y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];

plot(x,y,'*')

xlabel 'x轴'

ylabel 'y轴'

title '散点图'

.

hold on

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

end;

B=[0 0 0 0];

for j=1:n+1

for i=1:m+1

B(j)=B(j)+y(i)*x(i)^(j-1)

end

end

B=B';

a=inv(A)*B;

x=[-1.0:0.0001:2.0];

z=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3;

plot(x,z)

.

legend('离散点','y=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3') title('拟合图')