klobuchar2.asv

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

% function klobuchar(t,fm)
clear;
t=linspace(0,86000,100);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%变量初始值设置
%t=50000;
fm=linspace(-50,50,100);
%fm=linspace(20,80,71);
x0=[
5
10
0.07
0.01
0.005
0.000007
  
4000
200
15
2

2
0.5
7
-0.002073072

]

fx0=[
1
10
0.05
0.005
0.002
0.0000005
  
5000
100
20
1

1
-3
-0.5
-0.01073072

]
A1=x0(1,1);%%%%%%%%%%%%%更新x0的值
B=x0(2,1);
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);

 fA1=fx0(1,1);%%%%%%%%%%%%%更新x0的值
fB=fx0(2,1);
falpha0=fx0(3,1);
falpha1=fx0(4,1);
falpha2=fx0(5,1);
falpha3=fx0(6,1);
fbeta0=fx0(7,1);
fbeta1=fx0(8,1);
fbeta2=fx0(9,1);
fbeta3=fx0(10,1);  
fgama0=fx0(11,1);
fgama1=fx0(12,1);
fgama2=fx0(13,1);
fgama3=fx0(14,1);
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];
delta=zeros(1,100);
shigema=zeros(1,100);

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)-0.5)*2+f0(i,1);%%%%%%%%%%计算的初值
f(i,1)=It(fA1,fB,falpha0,falpha1,falpha2,falpha3,fbeta0,fbeta1,fbeta2,fbeta3,fgama0,fgama1,fgama2,fgama3,t(1,i),fm(1,i));%%%%%%%%%%计算的初值

end
for j=1:100;
    

x0=x0+dx;
    
    

H=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A1=x0(1,1);%%%%%%%%%%%%%更新x0的值
B=x0(2,1);
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);
x0=[A1;B;alpha0;alpha1;alpha2;alpha3;beta0;beta1;beta2;beta3;gama0;gama1;gama2;gama3];
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));
  x=2*pi*(t(1,i)-A3)/A4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   A1  
H(i,1)=1e-009; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  B
H(i,2)=-1e-009*fm(1,i);

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

if (A2>=0)&((abs(x)<1.57))

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   alpha1
if (A2>=0)&((abs(x)<1.57))

H(i,4)=fm(1,i)*cos(2*pi*(t(1,i)-A3)/A4);
else H(i,4)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    alpha2
if (A2>=0)&((abs(x)<1.57))
    
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)&((abs(x)<1.57))

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)&((abs(x)<1.57))&(A4<172800)

   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)&((abs(x)<1.57))&(A4<172800)


      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)&((abs(x)<1.57))&(A4<172800)

   % 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)&((abs(x)<1.57))&(A4<172800)

   % 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)&(A4/4>abs(t-A3))

    %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)&((abs(x)<1.57))


  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)&((abs(x)<1.57))

    %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)&((abs(x)<1.57))

    %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;
 

  delta(1,j)=sqrt((y-H*x0)'*(y-H*x0)/(1000-15));
if j>=2
  shigema(1,j)=abs(delta(1,j)-delta(1,j-1))/delta(1,j-1);
end
  if (j>=2)&((abs(delta(1,j)-delta(1,j-1))/delta(1,j-1))<0.01)

   %break;
end  
    

end

ff=f;