ぼんさい塾

ぼんさいノートと補遺に関する素材や注釈です.ミスが多いので初稿から1週間を経た重要な修正のみ最終更新日を残しています.

単振子の運動 (5)

2011-09-28 21:49:11 | 暮らし
phys.pdf
phys-s.pdf
記事一覧
 
                                  sinθの3次近似による計算

微分方程式

    d2x/dt2 + a x + b x2 + c x3 = 0

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

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


補足1: 2次のシンプレクティック法の 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)

補足2: 上図の 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


単振子の運動 (番外)

2011-09-26 21:41:18 | 暮らし

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)を変えたときの最初の5周期のピーク値(極大点近傍の3点を通る放物線の極値で近似)を次に示します.

    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 を半分にすると誤差はほぼ1桁減少しています.


単振子の運動 (4)

2011-09-24 19:26:42 | 暮らし
phys.pdf
phys-s.pdf
記事一覧

                          [2] 双線形化法による近似(3.4節)

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

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

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

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

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

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

補足3: 双線形化法はある種の非線形微分方程式の厳密解を求める強力な技法です.厳密解を探求する姿勢は差分化するときにも現れています.これに対して,[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)

2011-09-23 15:50:38 | 暮らし
phys.pdf
phys-s.pdf
記事一覧
 
                                  [2] 差分方程式の導出

後回しにしていた保存量

    H(t) = (1/2)(dx/dt)2 + (1/2) a x2 + (1/3) b x3 +(1/4) c x4

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

    c1 + c2 = a/4,   c3 = b/6,   c4 = c/6

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

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

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

補足3: 中点法では (1/2)(dx/dt)2 + (1/3) b x3 が保存量のときも x(nΔt + Δt/2)3 を (xn3 + xn+13) / 2 で近似するので,Δt → 0 の極限でしか保存量に一致しません.


単振子の運動 (2)

2011-09-21 15:36:05 | 暮らし

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

2011/09/26 本文・補足を修正

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

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

    d2x/dt2 + a x + b x2 + c x3 = 0

における保存量を

    H(t) = (1/2)(dx/dt)2 + (1/2) a x2 + (1/3) b x3 +(1/4) c x4

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

    dθ/dt = ∂H/∂ω = dEK/dω,    dω/dt = -∂H/∂θ = -dEU/dθ

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

(1) 1次のシンプレクティック法の場合,まず ωn+1 = ωn - Δt sin θn を求め,θn+1 = θn + Δt ωn+1 とします.
(2) 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 の場合で説明すると,a2 = 0, b2 = 1/2, a1 = 1, b1 = 1/2 です(A はθの更新,B はωの更新).係数 ai, bi は θn+1,ωn+1 のテイラー展開の誤差が O(tm+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
  一般化速度の変換則は接ベクトルの変換則と同じであり、一般化運動量の変換則は余接ベクトルの変換則と同じである。
  ※全然理解できません.