%用离散正交多项式求三次拟合多项式
% 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;
