1.軟體版本
matlab2013b
2.本算法理論知識
太陽黑子是人們最早發現也是人們最熟悉的一種太陽表面活動。因為太陽内部磁場發生變化，太陽黑子的數量并不是固定的，它會随着時間的變化而上下波動，每隔一定時間會達到一個最高點，這段時間就被稱之為一個太陽黑子周期。太陽黑子的活動呈現周期性變化是由施瓦貝首次發現的。沃爾夫 (R.Wolfer)繼而推算出11年的周期規律。實際上，太陽黑子的活動不僅呈11年的周期變化，還有海耳在研究太陽黑子磁場分布時發現的22年周期；格萊斯堡等人發現的80年周期以及蒙德極小期等。由于太陽黑子的活動規律極其複雜，時至今日科學家們仍在努力研究其内在的規律和特性。事實上，對太陽黑子活動規律的研究不僅具有理論意義，而且具有直接的應用需求。太陽黑子的活動呈現周期性變化的，沃爾夫(R.Wolfer)根據在過去的288 年(1700年~1987 年)間每年太陽黑子出現的數量和大小的觀測資料推算出11 年的周期規律。我們利用Matlab強大的資料處理與仿真功能，對Wolfer數進行功率譜密度分析進而可以得到對太陽黑子活動周期的結論。
使用的線性模型為
這個部分的MATLAB代碼為：
現象模型的預測結構如下：
将上面的圖放大之後，得到的局部估計效果如下圖所示：
然後通過多變量的最小二乘法進行模型預測。引入多個變量作為參數估計，這裡，将最小二乘法的估計函數為：
由于從曆年的資料可知，雖然總體上資料呈現出周期波形，但是對于每個周期中的波形，其中的局部圖中會出現部分高頻周期量，是以在估計的時候，這裡加上更高頻率的分量。使用和上面所講的最小二乘法進行上述參數的估計。
這個步驟，得到的仿真結果如下圖所示：
3.部分源碼
clc;
clear;
close all;
warning off;
load Sunspot.txt
%YEAR MON SSN DEV
%畫出導入的活動資料
YEAR = Sunspot(:,1);
MON = Sunspot(:,2);
SSN = Sunspot(:,3); %sun spot number
DEV = Sunspot(:,4); %标準偏差表
figure;
subplot(211);plot(SSN,'b.');
title('SSN');
subplot(212);plot(DEV,'b.');
title('DEV');
addpath 'funcs\'
%step1:Devise and/or employ methods to ?nd the best frequency for the sunspot cycle.
%step1:Devise and/or employ methods to ?nd the best frequency for the sunspot cycle.
%step1:Devise and/or employ methods to ?nd the best frequency for the sunspot cycle.
%找到一種較好的方法計算周期
K = 20;%前22年資料去除
Start = 12*K+1;
mainPeriod = func_cycle_cal(SSN(Start:end));
mainPeriod
%Step2:Predict the time of the next solar maximum and minimum.
%Step2:Predict the time of the next solar maximum and minimum.
%Step2:Predict the time of the next solar maximum and minimum.
%建構預測模型，預測後面的最大值和最小值對應的時間
%建構預測模型，預測後面的最大值和最小值對應的時間
c0 = [1 1 1 1 1 1 1 1]';%直接給出被辨識參數的初始值,即一個充分小的實向量
p0 = 10^6*eye(8,8); %直接給出初始狀态P0，即一個充分大的實數機關矩陣
t = 1/(12):1/(12):length(SSN)/(12);
yc = cos(2*pi*t/mainPeriod);
ys = sin(2*pi*t/mainPeriod);
yc2 = cos(2*pi*t/(2*mainPeriod));
ys2 = sin(2*pi*t/(2*mainPeriod));
yc3 = cos(2*pi*t/(3*mainPeriod));
ys3 = sin(2*pi*t/(3*mainPeriod));
for k=1:length(SSN); %開始求K
h1 =[1,k,yc(k),ys(k),yc2(k),ys2(k),yc3(k),ys3(k)]';
x = h1'*p0*h1 + 99;
x1 = inv(x); %開始求K(k)
k1 = p0*h1*x1; %求出K的值
d1 = SSN(k)-h1'*c0;
c1 = c0+k1*d1; %求被辨識參數c
e1 = c1-c0; %求參數目前值與上一次的值的內插補點
e2 = e1./c0; %求參數的相對變化
e(:,k) = e2; %把目前相對變化的列向量加入誤差矩陣的最後一列
c0 = c1; %新獲得的參數作為下一次遞推的舊參數
c(:,k) = c1; %把辨識參數c 列向量加入辨識參數矩陣的最後一列
p1 = p0-k1*k1'*[h1'*p0*h1+1]; %求出 p(k)的值
p0 = p1; %給下次用
end
%估計得到的太陽黑子活動曲線
for k=1:length(SSN); %開始求K
Y_predict2(k) = c(1,k) + c(2,k)*k + c(3,k)*yc(k) + c(4,k)*ys(k) + c(5,k)*yc2(k) + c(6,k)*ys2(k) + c(7,k)*yc3(k) + c(8,k)*ys3(k);
end
figure;
plot(Y_predict2,'b','LineWidth',2);hold on;
plot(SSN,'r');hold off;
legend('預測SSN','實際SSN');
Predict_Len = 100;%定義預測時間長度
t = 1/(12):1/(12):(length(SSN)+Predict_Len)/(12);
yc = cos(2*pi*t/mainPeriod);
ys = sin(2*pi*t/mainPeriod);
yc2 = cos(2*pi*t/(2*mainPeriod));
ys2 = sin(2*pi*t/(2*mainPeriod));
yc3 = cos(2*pi*t/(3*mainPeriod));
ys3 = sin(2*pi*t/(3*mainPeriod));
ind = 0;
for k=1:length(SSN)+Predict_Len; %開始求K
if k <= length(SSN)
Y_predict3(k) = c(1,k) + c(2,k)*k + c(3,k)*yc(k) + c(4,k)*ys(k) + c(5,k)*yc2(k) + c(6,k)*ys2(k) + c(7,k)*yc3(k) + c(8,k)*ys3(k);
else
%Y_predict3(k) = c(1,end) + c(2,end)*k + c(3,end)*yc(k) + c(4,end)*ys(k) + c(5,end)*yc2(k) + c(6,end)*ys2(k) + c(7,end)*yc3(k) + c(8,end)*ys3(k);
c0 = mean(c(1,end:end));
c1 = mean(c(2,end:end));
c2 = mean(c(3,end:end));
c3 = mean(c(4,end:end));
c4 = mean(c(5,end:end));
c5 = mean(c(6,end:end));
c6 = mean(c(7,end:end));
c7 = mean(c(8,end:end));
Y_predict3(k) = c0 + c1*k + c2*yc(k) + c3*ys(k) + c4*yc2(k) + c5*ys2(k) + c6*yc3(k) + c7*ys3(k);
ind = ind + 1;
Ys(ind) = Y_predict3(k);
end
end
figure;
plot(Y_predict3,'b','LineWidth',2);hold on;
plot(SSN,'r');hold off;
legend('預測SSN','實際SSN');
grid on;
%根據預測結果得到下次太陽黑子活動高峰和低峰的時間
%前一次高峰日期為
XX = 59;
[Vmax1,Imax1] = max(Ys);
[Vmax2,Imax2] = max(SSN(length(SSN)-XX:length(SSN)));%3100~3160
if Vmax1 > Vmax2
II = Imax1;
MM = Vmax1;
time = (length(SSN) + II-3019);%原資料的最後一個月份+預測後的最大值 - 前一個高峰日期
else
II = Imax2;
MM = Vmax2;
time = (length(SSN) + (XX-II)-3019);%原資料的最後一個月份+預測後的最大值 - 前一個高峰日期
end
Years=time/12;
fprintf('下次高峰期日期為：%d',round(2000 + Years));
fprintf('年\n\n');
fprintf('最大值為：%2.2f\n\n\n\n',MM);
%計算下一次低谷值
[Vmin,Imin] = min(Ys);
fprintf('下次低峰期日期為：%d',round(2012 + Imin/12));
fprintf('年\n\n');
fprintf('最小值為：%2.2f\n\n',Vmin);
A16-13