e0=[31.0615   40.4060   77.1353   21.1143]; %原始序列
[m,n]=size(e0);
a_l=exp(-2/(n+1))            %計(jì)算級(jí)比覆蓋范圍
a_h=exp(2/(n+1))
e0=sqrt(e0);                 %對(duì)原始數(shù)據(jù)作兩次方根變換   
e0=sqrt(e0);
for i=2:n
    q(i)=e0(i-1)/e0(i);      %計(jì)算實(shí)際級(jí)比
end
q
for i=1:n                     %計(jì)算一次累加序列
    e1(i)=sum(e0(1:i));       
end
for k=1:n-1                     
    z1(k)=1/2*e1(k)+1/2*e1(k+1);
end
Y=e0(2:n)';                    %得出Y矩陣
B=[];
B(:,1)=-z1';
B(:,2)=ones(1);                 %得出B矩陣
aU1=(inv(B'*B))*B'*Y;            %計(jì)算a,U 
a1=aU1(1);
U1=aU1(2);
for i=1:n 
    ee1(i)=(e0(1)-U1/a1)*exp(-a1*(i-1))+U1/a1; %解偏微分方程
end
for i=1:n-1
    ee0(i+1)=ee1(i+1)-ee1(i);           %一次累減，得出預(yù)測(cè)值
end
ee0=ee0.^4;                              %作四次方還原