太阳黑子周期分析.docVIP

. . 太阳黑子周期分析 1：计算太阳黑子周期 1）、选取历年的太阳黑子数据 本次作业选取的是1700—1999年的太阳黑子数据。将数据导入matlab中，并绘制太阳黑子数随年份变化的关系曲线。如图1所示。 程序如下: clear load sunspot.dat year =sunspot(:,1); sunspot =sunspot(:,2); plot(year(1:300),sunspot(1:300),b.-); xlabel (years); ylabel(sunspot data); title(1700—1999年太阳黑子数是随年份变化的关系曲线 ); grid on 图1、太阳黑子数随年份的变化曲线 2）：利用功率谱密度函数分析周期 1、对已经得到的Wolfer数进行FFT变换分析它的变化规律，并作功率与频率的关系图。 y=fft (sunspot (1:300)); y(1)=[]; n=length(y); power =abs(y(1:n/2)).^2; q=1/2; f= (1:n/2)/(n/2)*q; plot(f, power); xlabel(周期/年);title(周期图); 运行结果如图2所示。 图2、太阳黑子的功率谱 为了清楚起见，取功率和频率的前50个分量作它的周期图，程序如下: plot(f(1:50),power(1:50)); xlabel(频率); 运行结果如图3所示。 图3、功率和频率的前50个分量的周期图 2、确定太阳黑子的活动周期，画出功率与周期的关系图。程序如下： T=1./f; plot (T, power); axis ([0 50 0 7e+6]); %X轴范围是0-50，Y轴范围是0-7*10^6 xlabel (周期);ylabel(功率); grid on %在功率与周期的关系图上标出功率的最高点，该位置对应的周期即为太阳黑子活动的周期。程序如下： hold on index=find(power==max(power)); m=num2str(T(index)); plot(T(index),power(index),r.,MarkerSize,25); text(T(index)+2,power(index),[T=,m]); hold off 运行结果如图4所示： 图4、太阳黑子周期图 运用功率谱方法计算出太阳黑子的活动周期为T=11.0741，这与Wolfer得出的11年的周期规律基本一致，说明实验方法是正确的。 2、利用ARMA模型，预测未来某年的太阳黑子数 1）、建立AR模型 选用二阶自回归模型AR（2），方程为： （1） 采用最小二乘法对参数、进行估计： （2） 模型残差方差： （3） 计算参数程序如下： x=zeros(298,2); for i=2:1:299 x(i-1,1)=sunspot(i); end for k=1:1:298 x(k,2)=sunspot(k); end y=zeros(298,1); for t=3:1:300 y(t-2,1)=sunspot(t); end A=x; B=x*x; C=inv(B); D=C*A*y 运行得 D = 1.4867 -0.5981 即 带入公式（3）解得 求解程序如下： syms s m=0; s=sunspot (1:300); for i=3:300; m=m+(s(i)-1.4867*s(i-1)+0.5981*s(i-2))^2; end n=m/298 n = 364.1380 解得 故得到AR（2）模型方程是： 其中 2）、用上述AR（2）模型进行检验并预测 Sunspot(1998)=64.3, Sunspot(1999)=93.3, Sunspot(2000)=119.6 利用上述AR（2）模型计算得： Sunspot(2000)=1.4867*93.3-0.5981*64.3=100.2513 误差率=(119.6-100.2513）/119.6=16.18%; Sunspot(2004)=40.4, Sunspot(2005)=29.8, Sunspot(2006)=15.2 利用上述AR（2）模型计算得： Sunspot(2006)=1.4867*29.8-0.5981*40.4=20.1404; 误差率=（15.2-20.1404）/15.2=32.5% 经验证，AR（2）模型对之前所用数据的拟合程度是很好的，但是对后面年份的预测存在一定的