klobuchar.m

迭代实现非线性函数（电离层延迟模型）的参数拟合问题

% function klobuchar(t,fm)
t=linspace(100,70000,100);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%变量初始值设置
%t=50000;
fm=rand(1,100)*100;
%fm=linspace(20,80,71);
A1=5;
B=0.0004545454545;
alpha0=0.36;
alpha1=0.0441014;
alpha2=0.001;
alpha3=0.0003434534;
beta0=10000;
beta1=1000;
beta2=7;
beta3=3;
gama0=10;
gama1=0.6;
gama2=0.08;
gama3=0.01;

y=zeros(100,1);
H=zeros(100,14);
dx=zeros(14,1);
x00=zeros(14,100);
x0=[A1;B;alpha0;alpha1;alpha2;alpha3;beta0;beta1;beta2;beta3;gama0;gama1;gama2;gama3];
for i=1:100;
f0(i,1)=It(A1,B,alpha0,alpha1,alpha2,alpha3,beta0,beta1,beta2,beta3,gama0,gama1,gama2,gama3,t(1,i),fm(1,i));%%%%%%%%%%计算的初值
f(i,1)=rand(1)*10+f0(i,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设置观测值f(X)
%f(i,1)=0.5+f0(i,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设置观测值f(X)
end
%f=[25 24 23 22 22 23 23 24 23 22 19 15 12 11 12 15 21 28 36 42 48 52 54 55 56 57 59 61 65 70 76 85 95 106 118 129 140 149 156 161 164 165 165 163 160 155 149 142 135 131 129 131 136 143 150 156 158 157 153 148 142 138 136 135 133 130 123 113 100 87 77]
%f=f';

for j=1:100;
    


H=0;
A1=x0(1,1);%%%%%%%%%%%%%更新x0的值
B=x0(2,1);
%B=0.0004545454545;
alpha0=x0(3,1);
alpha1=x0(4,1);
alpha2=x0(5,1);
alpha3=x0(6,1);
beta0=x0(7,1);
beta1=x0(8,1);
beta2=x0(9,1);
beta3=x0(10,1);
gama0=x0(11,1);
gama1=x0(12,1);
gama2=x0(13,1);
gama3=x0(14,1);

for k=1:14;%%%%%%%%%%%%%%%%%%%储存每次迭代得到的x0
x00(k,j)=x0(k,1); 
end


for i=1:100;
    
 
A2=a2(alpha0,alpha1,alpha2,alpha3,fm(1,i))
A3=a3(gama0,gama1,gama2,gama3,fm(1,i))
A4=a4(beta0,beta1,beta2,beta3,fm(1,i))
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   A1  
H(i,1)=1; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  B
H(i,2)=t(1,i)-86400;
if  (0<=t(1,i))&(t(1,i)<=21600)
H(i,2)=t(1,i);
end
if  (21600<t(1,i))&(t(1,i)<72000)
    H(i,2)=t(1,i)-72000;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  alpha0

if (A2>=0)&(21600<t(1,i))&(t(1,i)<72000)

H(i,3)=cos(2*pi*(t(1,i)-A3)/A4);
else H(i,3)=0;
end

%if (21600<t<72000)
%H(i,3)=cos(2*pi*(t-A3)/A4);
%H(i,3)=cos(2*pi*(t-A3)/A4);
%else H(i,3)=0;
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   alpha1
if (A2>=0)&(21600<t(1,i))&(t(1,i)<72000)

%H(i,4)=fm(1,i)*cos(2*pi*(t-A3)/A4);
H(i,4)=fm(1,i)*cos(2*pi*(t(1,i)-A3)/A4);
else H(i,4)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    alpha2
if (A2>=0)&(21600<t(1,i))&(t(1,i)<72000)

%H(i,5)=fm(1,i)*fm(1,i)*cos(2*pi*(t-A3)/A4);
H(i,5)=fm(1,i)*fm(1,i)*cos(2*pi*(t(1,i)-A3)/A4);
else H(i,5)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     alpha3
if (A2>=0)&(21600<t(1,i))&(t(1,i)<72000)

%H(i,6)=fm(1,i)*fm(1,i)*fm(1,i)*cos(2*pi*(t-A3)/A4);
H(i,6)=fm(1,i)*fm(1,i)*fm(1,i)*cos(2*pi*(t(1,i)-A3)/A4);
else H(i,6)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    beta0
if (A4>=72000)&(21600<t(1,i))&(t(1,i)<72000)

   % H(i,7)=2*A2*pi*(t-A3)*sin(2*pi*(t-A3)/A4)/A4/A4;
   H(i,7)=2*A2*pi*(t(1,i)-A3)*sin(2*pi*(t(1,i)-A3)/A4)/A4/A4;
else H(i,7)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    beta1
if (A4>=72000)&(21600<t(1,i))&(t(1,i)<72000)

   % H(i,8)=2*A2*pi*(t-A3)*sin(2*pi*(t-A3)/A4)/A4/A4*fm(1,i);
      H(i,8)=2*A2*pi*(t(1,i)-A3)*sin(2*pi*(t(1,i)-A3)/A4)/A4/A4*fm(1,i);
else H(i,8)=0;
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    beta2
if (A4>=72000)&(21600<t(1,i))&(t(1,i)<72000)

   % H(i,9)=2*A2*pi*(t-A3)*sin(2*pi*(t-A3)/A4)/A4/A4*fm(1,i)*fm(1,i);
   H(i,9)=2*A2*pi*(t(1,i)-A3)*sin(2*pi*(t(1,i)-A3)/A4)/A4/A4*fm(1,i)*fm(1,i);
else H(i,9)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    beta3
if (A4>=72000)&(21600<t(1,i))&(t(1,i)<72000)

   % H(i,10)=2*A2*pi*(t-A3)*sin(2*pi*(t-A3)/A4)/A4/A4*fm(1,i)*fm(1,i)*fm(1,i);
   H(i,10)=2*A2*pi*(t(1,i)-A3)*sin(2*pi*(t(1,i)-A3)/A4)/A4/A4*fm(1,i)*fm(1,i)*fm(1,i);
else H(i,10)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  gama0 
if (43200<=A3)&(A3<=55800)&(21600<t(1,i))&(t(1,i)<72000)

    %H(i,11)=2*pi*A2*sin(2*pi*(t-A3)/A4);
     H(i,11)=2*pi*A2*sin(2*pi*(t(1,i)-A3)/A4);
else H(i,11)=0;
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  gama1
if (43200<=A3)&(A3<=55800)&(21600<t(1,i))&(t(1,i)<72000)

     %H(i,12)=2*pi*A2*sin(2*pi*(t-A3)/A4)*fm(1,i);
    H(i,12)=2*pi*A2*sin(2*pi*(t(1,i)-A3)/A4)*fm(1,i);
else H(i,12)=0;
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  gama2
if (43200<=A3)&(A3<=55800)&(21600<t(1,i))&(t(1,i)<72000)

    %H(i,13)=2*pi*A2*sin(2*pi*(t-A3)/A4)*fm(1,i)*fm(1,i);
    H(i,13)=2*pi*A2*sin(2*pi*(t(1,i)-A3)/A4)*fm(1,i)*fm(1,i);
else H(i,13)=0;
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  gama3
if (43200<=A3)&(A3<=55800)&(21600<t(1,i))&(t(1,i)<72000)

    %H(i,14)=2*pi*A2*sin(2*pi*(t-A3)/A4)*fm(1,i)*fm(1,i)*fm(1,i);
     H(i,14)=2*pi*A2*sin(2*pi*(t(1,i)-A3)/A4)*fm(1,i)*fm(1,i)*fm(1,i);
else H(i,14)=0;
end




end

H;
for i=1:100;
f0(i,1)=It(A1,B,alpha0,alpha1,alpha2,alpha3,beta0,beta1,beta2,beta3,gama0,gama1,gama2,gama3,t(1,i),fm(1,i));
%f(i,1)=rand(1)*100+f0(i,1);
end
y=f-f0;
dx=inv(H'*H)*H'*y;
x0=x0+dx;
end