単振子の運動 (5)

phys.pdf phys-s.pdf 記事一覧

sinθの３次近似による計算

微分方程式

    d²x/dt² + a x + b x2 + c x3 = 0

で x = θ, a = 1, b = 0, c = 1/6 とおいて d²θ/dt² = －sinθ を近似したときの計算式を上図に示します．また，２次のシンプレクティック法およびこの式のExcel による計算式を補足１，補足２に示します．これらは「単振子の運動 (1)」で示した式より簡単ですが，高次の近似であるため良好な結果が得られます．中点法では A1 を 0.3 にすると急速に減衰します．なお，中点法の「=D2+J2」「=E2+K2」等を「=D2」「=E2」等に変更すると２次の陽的ルンゲクッタ法による計算になります（このとき振幅は増大します）．

中点法（A1: 0.1） では  (B570, B571, B572)=(-0.993750228, -0.995159253, -0.988191219)
補足１ (A1: 0.3) では  (B189, B190, B191)=(-0.94491644,  -0.998367653, -0.976165956)
補足２ (A1: 0.3) では  (B190, B191, B192)=(-0.904983076, -0.992878482, -0.980457274)

補足１： ２次のシンプレクティック法の Excel による計算式．

A1: 0.1, B1: 1, C1: 0, A2: =A1/2

D101: =C101-A$2*SIN(B101)
B102: =B101+A$1*D101
C102: =D101-A$2*SIN(B102)

補足２： 上図の Excel による計算式．

A1: 0.1, A2: =A1*A1, A3: =1+A2/4, B1: 1, B2: =COS(A$1)

C102: =A$3+A$2*B102*B101/6
D102: =(2*B102-B101)-A$2*(2*B102+B101)/4
B103: =D102/C102

単振子の運動 (番外)

phys.pdf phys-s.pdf 記事一覧 9/30 全文差し替え

新B1: 0 旧B1 ＝ 新C1: 1 旧C1 ＝ 新D1: 0         新B2    : =B1+1         新B3以降: 上のセルのコピー

予告どおり 9/26 の記事を削除しました．代わりに Excel による計算について補足します．

「単振子の運動 (1)」で「Δt を変えても近似解があまり変わらなければほぼ収束しているとして扱えば気楽．」と述べましたが，Δt を半分にすると計算の行数は倍増します．しかし，例えば B201, C201 の値をテキスト化してB1，C1 に代入するだけで続きの計算を Excel がしてくれます．これに先立って，列 B の前に列を挿入して上図のように初期化しておけば，新B201，新C201, 新D201 の値を新B1，新C1，新D1 に代入したときにどの時間帯の計算をしているか分かり易くなります（グラフを挿入しておけば自動的に更新されます）．

参考までに，A1（＝Δt）を変えたときの最初の５周期のピーク値（極大点近傍の３点を通る放物線の極値で近似）を次に示します．

    A1=0.05        A1=0.1         A1=0.2         A1=0.3
  0.999935685    0.999484752    0.995841288    0.985809387
  0.999871369    0.9989695      0.991709667    0.971757659
  0.999807054    0.998454257    0.987558227    0.957670336
  0.999742738    0.99793903     0.983423741    0.943533985
  0.999678422    0.99742382     0.979277147    0.929569289

Δt を半分にすると誤差はほぼ１桁減少しています．

単振子の運動 (4)

phys.pdf phys-s.pdf 記事一覧

[2] 双線形化法による近似（3.4節）

最後に双線形化法（いわゆる広田の方法）を適用したときの計算について述べます．以下では，これまでの Δt でなく原著どおりのδを用います（Δ_t は演算子でΔt とは別物）．上図は計算過程における主要な式

(1) 双線形演算子 D_t の定義．
(2) x = g/f のときの dx/dt．
(3) 単振子の運動方程式．
(4) ゲージ不変性のある連立方程式（α＝0 にできる）．
(5) テイラー展開による遅延演算子の表現．
(6) 双線形差分演算子（中央差分）．
(7) 連立方程式のゲージ不変性を保持した差分化．
(8) 通常の中央差分が満たす式．
(9) 得られた差分方程式（a₁ + a₂ = a, b₁ + b₂ = b, Δ_t² x(t) = (x_n+1－2x_n＋x_n－1)/δ²）)．

を示していますが，原著をお持ちでない方には説明し切れません（説明するつもりで式番号を付けたのですが・・・）．とりあえず最終的に得られた式 (9) を見てください．青色で表示した部分は式を短くするために表現を変えています．つぎのことはすぐに分かります．

(a) 時間に関する対称性から x_n+1，x_n-1 のテイラー展開の１次の項が相殺される．
(b) 陰解法で解かねばならない．

補足１： (8) では x³ を x_n+1x_nx_n-1 で近似しています．中点法の (x_n+1 + x_n)³/8 による近似とは異なりますが，δが大きいと大域的な解軌道に厳密解とのずれが目立ってくることは中点法と同様であると思われます （pp.112-113 で保存量について説明されていますが，微分方程式の保存量との関係を理解できていません）． また x(t) の多項式でなく－sin x(t) のような式になると別途検討が必要になります．

補足２： 標本点での差分化を考えると奇数次の差分を標本値で表現するために相加平均をとる等の対策が必須です（※）．双線形化法では自由度を持たせていますが，中点法のように標本点の間での微分係数を差分化するだけでも簡単に時間軸での対称性を考慮した近似を実現できます．
※ (5) の Δ_tf(t) が標本点でない所の値を用いていることに注意． [7] の MP3 のような単純な差分化では Example 2 で例示した問題（spurious pole の発生）が生じます．

補足３： 双線形化法はある種の非線形微分方程式の厳密解を求める強力な技法です．厳密解を探求する姿勢は差分化するときにも現れています．これに対して，[1] に示した近似計算法は汎用性はありますが厳密解の検討には無力です．

[11] 広田の方法 - Wikipedia
  http://ja.wikipedia.org/wiki/%E5%BA%83%E7%94%B0%E3%81%AE%E6%96%B9%E6%B3%95
[12] 上野喜三雄，ソリトンが開く新しい数学，ISBN4-00-006504-1．
  http://www.iwanami.co.jp/.BOOKS/00/1/0065040.html
  「高校で学ぶ程度の・・・」ということですが難しい話が好きな人向き．この pp.64-65 に「広田の方法」が紹介されています．

単振子の運動 (3)

phys.pdf phys-s.pdf 記事一覧

[2] 差分方程式の導出

後回しにしていた保存量

    H(t) = (1/2)(dx/dt)² + (1/2) a x² + (1/3) b x³ +(1/4) c x⁴

の差分化を紹介します．ポイントは「保存量も時間反転に対して不変に」すべきであるとして，上図のように離散化します（記号はこのブログに合わせて変更しています）．c_i は

    c₁ + c₂ = a/4,   c₃ = b/6,   c₄ = c/6

となる定数です．上図の中括弧内の式を 0 とおく差分方程式の解は H_n+1 = H_n を満足します．

補足１： [2] p.89 に「失敗の原因は・・・振動方程式は・・・したがって差分方程式も・・・（ t → －t または＋δと－δの入れ替え）によって不変であるべきであり」と書かれていますが，微分可能な点では右極限（Δt → ＋0）と左極限（Δt → －0）が一致するので，δに関しては振動の場合に限りません．x_n+1 と x_n の混在についてはδだけでも説明できるような気がします．

補足２： 保存量が (1/2){ω(t)}² + {1 － cos θ(t)} のときに同様の技法 （因数分解できる形に修正） を適用することは難しそうです．なお p.91 の「・・・の２乗以上の項が現れ，解の一般性が失われるからである．」の意味は理解できていません．

補足３： 中点法では (1/2)(dx/dt)² + (1/3) b x³ が保存量のときも x(nΔt + Δt/2)³ を (x_n³ + x_n+1³) / 2 で近似するので，Δt → 0 の極限でしか保存量に一致しません．

単振子の運動 (2)

phys.pdf phys-s.pdf 記事一覧 2011/09/26 本文・補足を修正

[8] 指数作用素のフラクタル分割

広田先生の [2] には「ハミルトニアン」や「シンプレックティック」の語句が見当たりませんが非線形振動子の方程式

    d²x/dt² + a x + b x² + c x³ = 0

における保存量を

    H(t) = (1/2)(dx/dt)² + (1/2) a x² + (1/3) b x³ +(1/4) c x⁴

とおいて H(t) の差分化を説明しています．いうまでもなく「H」はハミルトニアンを意味していますが，一般化座標や一般化運動量の計算はありません．このため，[2] の紹介は次回にまわして，とりあえず [8] で教わったシンプレックティック数値積分の式を紹介します．運動方程式が dθ/dt = ω, dω/dt = －sin θのときの運動エネルギー E_K = (1/2)ω² はωだけの式，位置エネルギー E_U = (1－cosθ) はθだけの式です．シンプレクティック法では

    dθ/dt = ∂H/∂ω = dE_K/dω,    dω/dt = －∂H/∂θ = －dE_U/dθ

であることに着目してθ，ωを求めるのですが，時間発展のさせかたに特徴があります．

(1) １次のシンプレクティック法の場合，まず ω_n+1 = ω_n － Δt sin θ_n を求め，θ_n+1 = θ_n + Δt ω_n+1 とします．
(2) ２次のシンプレクティック法の場合は，ω_{n+ 1/2} = ω_n － (Δt/2) sin θ_n, θ_n+1 = θ_n + Δt ω_{n+ 1/2}, ω_n+1 = ω_{n+ 1/2} － (Δt/2) sin θ_n+1 です．
(3) m 次（m ≧ 3）の場合は上図に示されるような指数作用素のフラクタル分解を用います．最初の式の右辺をm = 2 の場合で説明すると，a₂ = 0, b₂ = 1/2, a₁ = 1, b₁ = 1/2 です（A はθの更新，B はωの更新)．係数 a_i, b_i は θ_n+1，ω_n+1 のテイラー展開の誤差が O(t^m+1) になるように定めます．

補足： f(t+Δt) を遅延作用素 Exp(Δt d/dt) = Σ_{k=0, ∞} (1/k!)(Δt d/dt)^k を用いて f(t+Δt) = Exp(Δt d/dt)f(t) と表わすことはよく行われますが，シンプレクティック法ではθだけを遅延させる作用素とωだけを遅延させる作用素が混在するため式が複雑になっています．

[8] Symplectic数値積分法
  http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/SYMP/index.html
  この文章は筆者が大学3年生の時に書いたものですので、・・・
  ※拍手!!
[9] 2 シンプレクティック法
  http://www.artcompsci.org/~makino/papers/ode/node3.html
  多変数だったら、もちろんまず速度の全要素を更新してから位置を更新しないといけないわけね。
  ※質量は定数ですから速度は運動量のようなものです．
[10] シンプレクティック幾何学 - Wikipedia
  http://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%B3%E3%83%97%E3%83%AC%E3%82%AF%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E5%B9%BE%E4%BD%95%E5%AD%A6
  一般化速度の変換則は接ベクトルの変換則と同じであり、一般化運動量の変換則は余接ベクトルの変換則と同じである。
  ※全然理解できません．