step_2.m

多项式拟和阶跃函数的小程序！对于多项式反演有一定的启发

clear all;close all;clc;
format long; %输出精度为1.0e+015
%原始数据

z = [50 48 46 44 42 40 38 36 34 32 30 29 28 27 26 25 24 23 22 21 20 ...
19 18 17 16 15 14 13 12 11 10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 ... 
3 2.5 2 1.5 1 0.5];
d = [28.5685 27.9427 26.3732 28.4689 28.0217 27.6379 28.6916 ... 
23.6709 24.3060 23.5071 25.7307 23.7017 24.3491 26.2110 ...
27.2651 25.1409 26.8494 24.2009 24.9812 23.9673 19.6244 ...
17.5687 21.4150 24.7871 24.9023 22.3715 21.9076 23.1678 ...
15.9585 4.6063 2.7828 2.3632 2.1839 3.4073 3.2262 2.6240 ...
2.3831 1.2161 1.5666 1.8210 1.9560 2.8652 1.7747 1.7908 ...
1.7763 1.7619 1.7474 1.7330 1.7185 1.6243];

%求解最小二乘解
N=50;X=0;X2=0;X3=0;X4=0;X5=0;X6=0;X7=0;X8=0;X9=0;X10=0;X11=0;X12=0;
X13=0;X14=0;X15=0;X16=0;X17=0;X18=0;X19=0;X20=0;X21=0;X22=0;
Di=0;  %GtG矩阵中的各元素(N,X,X1,...,X16)
%对Gtd矩阵中的列向量赋初值
M1=0;M2=0;M3=0;M4=0;M5=0;M6=0;M7=0;M8=0;M9=0;M10=0;M11=0;
for i=1:50;
    X=X+z(i);        %求z(i)从1到50的和
    X2=X2+z(i).^2;   %求z的平方和，从1到50，以下依此类推
    X3=X3+z(i).^3;
    X4=X4+z(i).^4;
    X5=X5+z(i).^5;
    X6=X6+z(i).^6;
    X7=X7+z(i).^7;
    X8=X8+z(i).^8;
    X9=X9+z(i).^9;
    X10=X10+z(i).^10;
    X11=X11+z(i).^11;
    X12=X12+z(i).^12;
    X13=X13+z(i).^13;
    X14=X14+z(i).^14;
    X15=X15+z(i).^15;
    X16=X16+z(i).^16;
    X17=X17+z(i).^17;
    X18=X18+z(i).^18;
    X19=X19+z(i).^19;
    X20=X20+z(i).^20;
    X21=X21+z(i).^21;
    X22=X22+z(i).^22;
    Di=Di+d(i);        %求观测数据的和
    %求Gtd矩阵中的9个列向量
    M=Di;M1=M1+z(i).*d(i);M2=M2+z(i).^2.*d(i);M3=M3+z(i).^3.*d(i);
    M4=M4+z(i).^4.*d(i);M5=M5+z(i).^5.*d(i);M6=M6+z(i).^6.*d(i);
    M7=M7+z(i).^7.*d(i);M8=M8+z(i).^8.*d(i);M9=M9+z(i).^9.*d(i);
    M10=M10+z(i).^10.*d(i);M11=M11+z(i).^11.*d(i);
end
GtG=[N X X2 X3 X4 X5 X6 X7 X8 X9 X10 X11;
     X X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12;
     X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13;
     X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14;
     X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15;
     X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16;
     X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17;
     X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18;
     X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19;
     X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20;
     X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21;
     X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21 X22];
A=inv(GtG);%求GtG的逆矩阵
Gtd=[M;M1;M2;M3;M4;M5;M6;M7;M8;M9;M10;M11];
m=A*Gtd    %求出模型
plot(z,d,'.')%将原始数据以'.'的形式在图形中显示
hold on
y1=m(1)+m(2).*z+m(3).*z.^2+m(4).*z.^3+m(5)*z.^4+m(6)*z.^5+m(7)*z.^6 ...
+m(8)*z.^7+m(9)*z.^8+m(10)*z.^9+m(11)*z.^10+m(12)*z.^11;
plot(z,y1,'r') %将反演出的模型和已知数据在同一图中显示，标为红色
%误差评价，相关系数:Rxy,拟合误差
%相关系数r
aver_y1=sum(y1)/50;
aver_d=sum(d)/50;
Lyy=sum((y1-aver_y1).^2);
Lxx=sum((d-aver_d).^2);
Lxy=sum((y1-aver_y1).*(d-aver_d));
r=Lxy/sqrt(Lxx*Lyy)
%拟合误差
err=sum((d-y1).^2)/50