离散正交的二次多项式matlab实现,月光软件站 - 编程文档 - 其他语言 - 用离散正交多项式求三次拟合多项式[MATLAB版本]...

%用离散正交多项式求三次拟合多项式

% x,y--表示原始数据的节点坐标

% w--表示权重系数

% N--表示要拟合的离散正交多项式的最高次数

% polyapproximate()--是自定义函数，可以求解多项式的系数

% 其返回值c为多项式系数，error为偏差平方和

x=[0 0.2 0.6 1.0 1.3 1.6 1.7 1.8 1.9 2.2 2.3 2.5 2.6 2.9 3.1 3.4 3.8 4.1 4.4 4.7 4.8 4.9 5.0 5.1 5.3];

nn=length(x);

for i=1:nn

w(i)=1;

end

y=[0 -2.5 -4.0 -5.7 -3.5 -2.0 -1.0 2.0 3.5 4.0 7.0 7.5 9.9 10.9 11.9 13.5 13.0 11.9 9.0 6.5 4.0 1.5 0.0 -2.5 -5.0];

N=3;%此处可取3 or 4.

[c,error]=polyapproximate(x,y,w,N)

t=0:0.1:5.3;

u=polyval(c,t);

plot(t,u,x,y,'+')

%自定义函数polyapproximate()，用来做离散正交多项式拟合

% 此函数的作用是做不同次数的离散正交多项式的拟合

% X，Y 为原始数据的坐标值矩阵

% w 为权重系数

% N 为离散正交多项式的最高次数

function [C,E]=polyapproximate(X,Y,w,N)

M=length(X);

for i=1:N+1

for j=1:i

if j~=i

P(i,j)=0;

else

P(i,j)=1;

end

end

end

S=0;

d(1)=0;

for i=1:M

d(1)=d(1)+w(i);

S=S+w(i)*X(i);

end

AF(1)=S/d(1);

P(2,1)=-AF(1);

for i=1:M

PX(i,1)=1;

PX(i,2)=X(i)-AF(1);

end

BA(1)=0;

for k=2:N+1

S=0;

dd=0;

for i=1:M

S=S+w(i)*X(i)*PX(i,k)*PX(i,k);

dd=dd+w(i)*PX(i,k)*PX(i,k);

end

d(k)=dd;

AF(k)=S/d(k);

BA(k-1)=d(k)/d(k-1);

P(k+1,1)=-AF(k-1)*P(k,1)-BA(k-1)*P(k-1,1);

for i=1:k-1

j=k-i+1;

if j>=k

t=0;

else

t=P(k-1,j);

end

P(k+1,j)=P(k,j-1)-AF(k-1)*P(k,j)-BA(k-1)*t;

end

for i=1:M

PX(i,k+1)=PX(i,k)*(X(i)-AF(k-1))-BA(k-1)*PX(i,k-1);

end

end

d(N+1)=0;

for i=1:M

d(N+1)=d(N+1)+w(i)*PX(i,N+1)*PX(i,N+1);

end

for i=1:N+1

FM=0;

for k=1:M

FM=FM+w(k)*Y(k)*PX(k,i);

end

gp(i)=FM/d(i);

end

for i=1:N+1

C(i)=0;

for j=i:N+1

C(i)=C(i)+gp(j)*P(j,i);

end

end

C=flipud(C');

%C=C'

U=0;

for i=1:M

U=U+w(i)*Y(i)*Y(i);

end

V=0;

for k=1:N+1

V=V+gp(k)*gp(k)*d(k);

end

E=U-V;