カオス

時間変化が非周期的で予想不可能な非線形な系の運動をカオス (Chaos) という．振り子などの機械的に振動する物体、対流する液体、気体、レーザーキャビティー、化学反応などいろいろな系にみられる現象である．カオスの本質は過去におこなった運動を再びおこさないことにある．カオス系の運動を記述する運動方程式ははっきりと未来の運動を予想する (deterministic equations) のであるが、実際の運動は全く規則性を欠き初期条件がほんの少し違うと時間が少したっただけでも全く違う軌道を描く．
　わかりやすく説明するために、いまある系がほんの少し違う初期条件で運動を始めたとする．もし、系がカオス的でないとすると、その初期における違いは時間と共に線形に増大する．カオス系はその違いが指数関数的に増大するのである．初期条件がちょと違えばもうその後の運動はあっと言う間に全く違ってしまう．この意味において、カオスでは予想不可能であるといっているのである．このことについて最初に言及したのは HenriPoincare (1854-1912) である．
    教科書と同じ方程式を数値的に解くことにする．ただし、全く同じに行ってもおもしろくないので、ここではNDSolve を使って解くことにする．
　長さ L の糸の先に質量 m の振り子がぶら下がり、速度に比例する摩擦力 
γ m L ((dθ)/dt)、それに外力 (Subscript[F, 0]) cos ((Subscript[ω, D]) t) を受けていると考える．ここで γ は定数で減衰パラメターと呼ぶ．そのときのニュートン の運動方程式は、
m L ((d^2)/(dt^2))θ + γ m L((dθ)/dt)+ W sinθ =  (Subscript[F, 0]) cos ((Subscript[ω, D]) t) 
とかける．W は重量、 (Subscript[ω, D]) は外力の角振動数である． ここで、系のパラメターの数を減らすために、振り子の自然振動数 (Subscript[ω , 0])(= (Sqrt[(W/(m L))]))   が１であると仮定する．時間はしたがって無単位とみなす．この仮定のもとで、運動方程式は次のようにかきかえられる．
ω'[t] + 1/q  ω[t] + Sin[θ[t]] == g Cos[ωD t]
ここで、q はクオリティーパラメター (quality parameter)、 g は外力振幅 (forcing amplitude)と呼ぶことにする． (g を重力の加速度と混同しないように．)

初期条件

正弦関数的な外力が振り子（減衰振動系）にかかっているとする．詳しくは教科書を参照．初期条件を決めて、運動方程式を解こう．例えば、

{θ0, ω0, q, g,ωD}={-.8,.7,4,1.5,Pi/6};

solution=NDSolve[{ω'[t]+1/q  ω[t]+Sin[θ[t]]==g Cos[ωD t], ω[t]==θ'[t], θ[0]==θ0, ω[0]==ω0}, { θ[t], ω[t]}, {t, 0, 100}]

ParametricPlot[Evaluate[{ θ[t], ω[t]}/.solution[[1]]], {t, 0,89},AxesLabel->{"θ[t]", "ω[t]"}]

となり、相空間における軌跡を描くことができる．ここで、相空間での運動の軌跡をさらに長時間にわたって描いてみよう．ただ角度は２Piの周期を持っているので、教科書のように（117ページ)theta の範囲は-PiとPiとの間で描いてみよう．

solution=NDSolve[{ω'[t]+1/q  ω[t]+Sin[θ[t]]==g Cos[ωD t], ω[t]==θ'[t], θ[0]==θ0, ω[0]==ω0}, { θ[t], ω[t]}, {t, 0, 420}, MaxSteps->10000]

ParametricPlot[Evaluate[{Mod[θ[t]+π,2π]-π,ω[t]}/.solution[[1]]],{t,0,420},AxesLabel->{"θ[t]","ω[t]"}]

これでは、詳細がわかりにくいのでListPlotを使って描き直そう．

ListPlot[Table[Evaluate[{Mod[θ[t]+π,2π]-π,ω[t]}/.solution[[1]]],{t,0,420,.03}],AxesLabel->{"θ[t]","ω[t]"}, PlotStyle->{PointSize[.001], Hue[1]}]

（教科書の口絵15を参照）．

さて、初期値を若干変えて、軌跡がどのように変化するか見るために、いろいろプロットを行ってみよう．

Do[{θ0, ω0, q, g,ωD}={-.8,.6+n/5,4,1.5,Pi/6}; 
solution=NDSolve[{ω'[t]+1/q  ω[t]+Sin[θ[t]]==g Cos[ωD t], ω[t]==θ'[t], θ[0]==θ0, ω[0]==ω0}, { θ[t], ω[t]}, {t, 0,24}];
graph[n]=ParametricPlot[Evaluate[{ θ[t], ω[t]}/.solution[[1]]], {t, 0,24},AxesLabel->{"θ[t]", "ω[t]"}, PlotStyle->Hue[.1n]], {n, 10}]

Manipulate[Show[graph[n]],{n,1,10,1}]

これらのグラフを一つにして描くと、

Show[graph[1],graph[2],graph[3],graph[4],graph[5],graph[6],graph[7],graph[8],graph[9],graph[10], PlotRange->All]

すなわち、初期値が少し違うと運動の軌跡は時間がたつと（上の例では外力の周期の２倍の時間までプロットした．）全く違うものになることが推察できる．

時間が０、半周期、１周期、１．５周期、２周期後の”位置”を描いてみよう．

datalist={};
{ q, g,ωD}={4,1.5,Pi/6};Do[Do[{θ0, ω0}={-1+n/20,-1+ m/20}; 
solution=NDSolve[{ω'[t]+1/q  ω[t]+Sin[θ[t]]==g Cos[ωD t], ω[t]==θ'[t], θ[0]==θ0, ω[0]==ω0}, { θ[t], ω[t]}, {t, 0,24}];
datalist=Prepend[datalist, Evaluate[Table[{ θ[t], ω[t]}/.solution[[1]], {t, 0, 24, 6}]]], {n, 20}], {m, 40}]

Manipulate[graph[n]=ListPlot[Evaluate[Table[datalist[[m,n]], {m,Length[datalist]}]], PlotStyle->Hue[.2n]], {n, 1,5,1}]

これらを一つにしてプロットしてみると、

Show[graph[1], graph[2],graph[3], graph[4],graph[5]]

それぞれのデータ点は相空間上に置いて、時間が０、半周期、１周期、１．５周期、２周期後の”位置”を表す点に対応する（上図は教科書の口絵１４に対応する）．初期値が少し違えば大きく拡散することと、ある集合上に（奇妙なアトラクターという）点が集まるであろうことを彷彿させる．

問1

上の例を参考に初期値を

{ q, g, ωD} = {1., 1., Pi/6};

や

{ q, g, ωD} = {2, 2.5, Pi/6};

として、同じ計算をし、プロットしてみなさい．結果の考察を行いわかりやすく説明して下さい．（教科書112ページからの記述を参照）．

ポアンカレ断面図

上の例では周期はt=12 である．位相がPi/3, 2Pi/3, Pi, 4/3 Pi, 5/3 Pi, 2Piの時の位相空間での粒子の位置を求めて試しにプロットしてみよう．教科書にもあるようにこれがポアンカレ断面図の例であり、入れ子構造をもつフラクタル図を得ることが出来る．

datalist1={};
{q,g,ωD}={4,1.5,(π/6)};
{θ0,ω0}={-.7,-1};
solution=NDSolve[{ω'[t]+((ω[t])/q)+Sin[θ[t]]==g Cos[ωD t],ω[t]==θ'[t],θ[0]==θ0,ω[0]==ω0},{θ[t],ω[t]},{t,0,12}];
Timing[Do[{θ0,ω0}={θ[t],ω[t]}/.solution[[1]]/.t->12;solution=NDSolve[{ω'[t]+((ω[t])/q)+Sin[θ[t]]==g Cos[ωD t],ω[t]==θ'[t],θ[0]==θ0,ω[0]==ω0},{θ[t],ω[t]},{t,0,12}];datalist1=Prepend[datalist1,Evaluate[Table[{Mod[θ[t]+π,2 π]-π,ω[t]}/.solution[[1]],{t,2,12,2}]]],{n,15000}]]

Manipulate[graph2[n] =ListPlot[Take[Table[datalist1[[m,n]], {m, Length[datalist1]}],15000],PlotStyle->{PointSize[.003], Hue[.16n]}], {n,1,6,1}]

Show[graph2[1],graph2[2],graph2[3],graph2[4],graph2[5],graph2[6]]

問２

上の図で位相とポアンカレ図とはどのような関係にあると推測できるか．

問３

上の計算で、時間をもっと長くとってもっとおおくの点をプロットとすると、図はどのように変化するか．

問４

上の計算で、初期値を

{θ0,ω0}={-.7,-1};

としたが、他の値を使うと結果はどうなるか．

問５

上の計算で、q, g の値を変えて計算をすると結果はどのようになるであろうか．（教科書を参照）．

問 ６

以上、上の考察で得られた結果をもとに、カオスに関する自分の意見をまとめよ。感想でもよい。