演習例題

前節のとおり、データの共分散行列を求め、それに関わる固有値問題を解くことでEOFsを求められることが分かった ここでは、例題に挑戦してみよう。

次に示す、三角関数の線形結合で糞れるデータをEOFsに分けて、第1と第2モードを錘ｦしよう。 データは時間$t$ を0から5$\pi$ の間、空間では0から$\pi$ の間で取得されるとする。次のような 4つの定在波を聞三角関数を想定し、その線形結合を議論する。

$$a_1=sin\left[\omega_1 t-k_1(x-\pi/2))+sin\left[\omega_1 t+k_1(x-\pi/2))+0.1R\right]$$ (16)

$$a_2=-0.5sin\left[\omega_2 t-k_2(x-\pi/4))-0.5sin\left[\omega_2 t+k_2(x-\pi/4))+0.05R\right]$$ (17)

$$a_3=0.25sin\left[\omega_3 t-k_3(x-\pi/6))+0.25sin\left[\omega_3 t+k_3(x-\pi/6))+0.025R\right]$$ (18)

$$a_4=-0.1sin\left[\omega_4 t-k_4(x-\pi/8))-0.1sin\left[\omega_4 t+k_4(x-\pi/8))+0.0125R\right]$$ (19)

ここで$k$ は波数、$\omega$ は周波数を聞。この4つのデータをMatlabで生成し其々、及び、それらの線形結合の変動を錘ｦしてみよう。

x=0:0.02:pi;
t=0:0.2:5*pi;
[t,x]=meshgrid(t,x);
k=[1 2 3 4];
omega(1)=1;
omega(2)=2;
omega(3)=3;
omega(4)=4;
a(:,:,1)=sin(omega(1).*t-k(1).*(x-pi/2))+sin(omega(1).*t+k(1).*(x-pi/2))+0.1.*rand(size(x));
a(:,:,2)=-0.5.*sin(omega(2).*(t-pi/4)-k(2).*(x-pi/4))-0.5.*sin(omega(2).*(t-pi/4)+k(2).*(x-pi/4))+0.05.*rand(size(x));
a(:,:,3)=0.25.*sin(omega(3).*(t-pi/6)-k(3).*(x-pi/6))+0.25.*sin(omega(3).*(t-pi/6)+k(3).*(x-pi/6))+0.025.*rand(size(x));
a(:,:,4)=-0.1.*sin(omega(4).*(t-pi/8)-k(4).*(x-pi/8))-0.1.*sin(omega(4).*(t-pi/8)+k(4).*(x-pi/8))+0.0125.*rand(size(x));
D=sum(a,3);

figure
for i=1:1:4
  subplot(2,2,i)
  pcolor(t,x,a(:,:,i))
  shading flat
  colorbar
  title(['a_' num2str(i)],'fontsize',15)
  xlabel('time','fontsize',15)
  ylabel('y','fontsize',15)
end
saveas(gcf,'waves4','png')
print(gcf,'-depsc','waves4')

figure
pcolor(t,x,D);shading flat;colorbar;
title('D','fontsize',15)
xlabel('time','fontsize',15)
ylabel('y','fontsize',15)
saveas(gcf,'sumwaves4','png')

このプログラムを名前をつけて保存して実効してみよう。 以下の様な翠謔ｩける。

Figure: 式(16-19)で生成される定在波的なデータを時間と空間の関数として各々示す。

Figure: Fig.(1)を線形結合したものを示す。

次に、データ $\mathbf{D}$ の共分散行列 $\mathbf{C}$ を求め、 $\mathbf{C}$ の 固有値 $\mathbf{L}$ と固有ベクトル $\mathbf{E}$ を求めよう。固有値問題はMatlabが常備する行列の 計算ツールの内、固有値問題の解法を実行する$eig$ を用いる。

C=D*D'./size(D,2);　% D'はDの転置行列 D*D'./Nは共分散行列
[E,L]=eig(C);

結果として得られた、固有値 $\mathbf{L}$ と元のデータ $\mathbf{D}$ を用いて、 各モードの振幅$a_i(t)$ を求める。具体的には式(15)を実際に計算する。 その後、得られた $\mathbf{A}=a_i(t)$ と固有ベクトル $\mathbf{E}$ から第1モード と第2モードの時空間困を抽出する。抽出は、元の式(6)を用いれば良い。

A=E'*D;
reD1=E(:,end)*A(end,:);
reD2=E(:,end-1)*A(end-1,:);

これで、reD1に第1モード、reD2に第2モードを格納した。 これらを次の様に錘ｦすると以下の様になる。

figure
pcolor(t,x,reD1)
shading flat
colorbar
title('1^{st} mode','fontsize',15)
xlabel('time','fontsize',15)
ylabel('y','fontsize',15)
saveas(gcf,'mode1st','png')
print(gcf,'-depsc','mode1st')

figure
pcolor(t,x,reD2)
shading flat
colorbar
title('2^{nd} mode','fontsize',15)
xlabel('time','fontsize',15)
ylabel('y','fontsize',15)
saveas(gcf,'mode2st','png')
print(gcf,'-depsc','mode2nd')

figure
pcolor(t,x,reD1+reD2)
shading flat
colorbar
title('sum of 1^{st} and 2^{nd} mode','fontsize',15)
xlabel('time','fontsize',15)
ylabel('y','fontsize',15)
saveas(gcf,'summode12','png')
print(gcf,'-depsc','summode12')

Figure: 抽出されたEOF第1モードの時空間困を示す。

Figure: 抽出されたEOF第2モードの時空間困を示す。

抽出された第1、2モードの線形結合は、元のデータ $\mathbf{D}$ (Fig. )に酷似している。これは最初の2つのモードが時空間変動に 支配的であることを示す。しかし若干元データと異なるのは、他のモードからの寄与が無視できない空である。 しかしながら、 EOFsは経験的な統計手法でしかなく、実際の機昨反映し得ないため、EOFsの解釈には常に注意が必要である。

Figure: 抽出されたEOF第1,2モード総和の時空間困を示す。