時系列分析

　２０１２年という映画が今年2009年の11月に公開される。私は、ノストラダムス世代なので、子供のころ、きっと１９９９年には世界は滅ぶと思っていたので、こういう映画には少し惹かれる。
なんでも、マヤ文明の暦が２０１２年の冬至に終わりを告げることになっていて、この日に人類が滅亡するという内容らしい。
　以前、マヤ文明は高度な暦をもっていて、旱魃などを正確に予言していたということを聞いたので少し調べてみた。しかし、マヤ文明は確かに高度な暦をもっていたらしいが、2012年の冬至に滅ぶという科学的根拠は全くないらしい。マヤ文明とあまり関係のない、太陽系がフォトンベルトに突入して、太陽が異常活動して滅ぶという荒唐無稽な説と合体して人類滅亡説ができあがっているようだ。
　２０１２年には日本の東京でも金環日食がみられるらしいので楽しみな年だが、こうした珍しい天体現象も予言と結び付けられているようだ。
　
世界が滅亡説は、隕石の落下、異常気象、メタンハイブレードの爆発、スーパ火山の噴火、太陽フレアなど様々なことが想定されて、映画になったりしているが、なかなか興味深い。

　それに関連して、地磁気が周期的に反転しているということを聞いたことがあり、人類が滅ぶとしたら、地磁気反転で、太陽風によって人類が被爆して死ぬということもありうるなあと考えたことがあったので、今一度、地磁気反転について少し調べてみた。
　地磁気は、ダイナモ理論により説明されるというのは、聞いていたが、東大名誉教授だった力武常次先生が立てられた力武モデルというのがすごい。
　これは、力武先生が、1957年発表した結合円盤ダイナモモデルというもので、2つの円盤型発電機があって、それぞれの発電機の円盤からブラシにより電流を取出し、コイルに流して、円盤の軸に電流を戻す電気回路で発電機を結合させることで、地磁気の反転を説明しようというものだ。実際、人工衛星の観測などで、通常の向きとは逆の磁束斑が認められるようになっており、この説は広く認められてきているようだ。

このモデルの電気回路の式は、

L1 dI1/dt + R1 I1 = ω1 M1 I2
L2 dI2/dt + R2 I2 = ω2 M2 I1
T1 dω1/dt = G1 - T1 M1 I1 I2
T2 dω2/dt = G2 - T2 M2 I1 I2

であらわされる。
ここで、L1,L2は自己インダクタンス、R1R2は電気抵抗、M１,M2は相互インダクタンス、T1,T2は円盤の回転の慣性能率、G1,G2は発電機を回すトルクである。

簡単のためL1=L2,R1=R2、M1=M2　、T1=T2,G1＝G2とおいて、式を書き直すと

L dI1/dt + R I1 = ω1 M I2
L dI2/dt + R I2 = ω2 M I1
T dω1/dt = G - M I1 I2
T dω2/dt = G - M I1 I2

ここで無次元化するために

I1 = J1 √（G / M）　, I2 = J2 √（G / M）
ω1 = Ω1 √（G L/T M) ,ω2 = Ω2 √（G L / T M)

t = τ √（T L / G M)

とすると

dJ1/dτ　+　R√（T/GML）J1　=　Ω1　J2
dJ2/dτ　+　R√（T/GML）J2　=　Ω2　J1
dΩ1/dτ = 1 -J1 J2
dΩ2/dτ = 1 - J1 J2

に変換され、
ここでμ=R√（T/GML）とおくと

dJ1/dτ + μ J1 = Ω1 J2
dJ2/dτ + μ J2 = Ω2 J1
dΩ1/dτ = 1 - J1 J2
dΩ2/dτ = 1 - J1 J2

J1=X,J2=Y,Ω1=Z　Ω1-Ω2=α(一定）とすると

dX/dτ = -μX + ZY
dY/dτ = -μY － αX + XZ
dZ/dτ = 1 - XY

に変換できる。
ようするにαとμを定数とした常微分方程式に変換される。
これは、カオスの教科書に出てくるあのローレンツ方程式に似ているが、実際カオス的な振る舞いをする。
ちなみにこのモデルはローレンツが決定論的非周期な流れとしてローレンツ方程式を示した1963年より6年も前に提唱されている。
恐るべし、力武常次教授。

R言語でプロットしてみると
イメージ 2

プログラムは以下のとおり

library(scatterplot3d)
mu=2
k=3
a=mu*(k*k-1/(k*k))
dt=0.01
x[1]<-0.09
y[1]<--0.3
z[1]<-18
for(i in seq(5000)){
	dx = (- mu*x[i] + z[i]*y[i]) * dt
	dy = (- mu*y[i] - a*x[i] + z[i]*x[i]) * dt
	dz = ( 1- x[i]*y[i] ) * dt
	x[i+1]=x[i]+dx
	y[i+1]=y[i]+dy
	z[i+1]=z[i]+dz
		
}


scatterplot3d(x, y, z, highlight.3d=TRUE, 
              col.axis="blue", col.grid="lightblue",
              main="rikitake model", type="l")

実際の地球はもっと多くの発電機が連結されたものとなるのだろうが、これにより、地磁気が互いの電流と磁場の影響を受けながら、反転していくことが説明できる。
図にｘ,ｙ,ｚの変化を示す。ｘ,yが反転しているのが見て取れる。
イメージ 3

par(mfrow=c(3,1))
plot(x,type="l")
abline(h=0)

plot(y,type="l")
abline(h=0)
plot(z,type="l")

カオス理論によれば、いつ地磁気反転がおきるかは周期的ではなく、いつ地磁気が消えてしまってもおかしくなく結構怖いが、短期予測は可能である。
地球の磁場解析はホットな研究テーマらしく、今も盛んに研究されている。
地磁気の反転する時期も予測できるようになるかもしれない。

なんかくだらないことを考えて雑学をためこんでしまった。

# 作図準備 x <- seq(0, 5, by=0.1) y <- seq(0, 1, by=0.02) plot( x, y, ylab="", type='n', main = 'Logistic Map - Bifurcation Diagram -') # 初期値 x0=0.1 #aを1から400まで0.01づつ変化させる for (j in 100:400){ a=j/100 x=x0 for(i in 1:100){ y=a*(1-x)*x if(i >=50){ #Xn n>=50の点を50点プロットする points(a,y,pch='.') } x=y } }