単振子の運動 (1)

phys.pdf phys-s.pdf 記事一覧 9/23 説明追加

[6]中点法の近似解法

まず原始的な中点法から復習します．簡単のため ω(nΔt) を ω_n と略記します．

   θ'(t) = ω(t),    ω'(t) = －sin θ(t)

とすると

   ω_n+1－ω_n+1＝－∫_{[nΔt, (n+1)Δt]} sin θ(t) dt

であり，右辺は中間値の定理で －sin θ(τ)Δt の形で表せます．このτを (n + 1/2)Δt で近似し，θ(τ) を (θ_n＋θ_n+1)/2 で近似するのが最も基本的な中点法で，上図の差分方程式は次のようになります．

   (θ_n+1－θ_n)/Δt = (ω_n＋ω_n+1)/2,    (ω_n+1－ω_n)/Δt = －sin((θ_n＋θ_n+1)/2)

このような中点法は陰解法なので，方程式が非線形の場合厳密解を求めることは一般に困難ですが，[6] の方法によって近似解を求めることができます(上図参照)．しかし，中点法自体が近似であり，Δt → 0 の極限で微分方程式の解に一致するというだけです．それでも「中点法を使えばよい」というのは次の理由によります：(1)微分は１次近似の極限であり，極限をとる前でも基本的な特徴（例えば解軌道の曲がり方）が現れていることが多い．(2)Δt を変えても近似解があまり変わらなければほぼ収束しているとして扱えば気楽．

補足１： 微分方程式が線形の場合，ラプラス変換のｓ平面の原点近傍が中点法で離散近似したｚ変換の z =１ の近傍に比較的精度よく対応します．インパルス応答の標本値を保存する変換（標準ｚ変換）があるにも拘わらず，そうでない変換（双１次変換）をフィルタの設計に用いるのはこのためです --- 虚軸上に極を持つような回路はフィルタに使えません．減衰振動の仕方も好ましく近似されていることが重要です．

補足２： sin θ(τ) を sin θ(nΔt + Δt/2) で近似するとエネルギーは保存されません（ほぼ保存するように強引に補正することは可能[1]）．[6] の方法による Excel の式を次に示します．

A1: 0.1, B1: 1, C1: 0, A2: =A1/2, A3: =A1*A1/2, A4: =A1*A1/4

B102: =D101+J101 
C102: =E101+K101
D102: =B102+A$1*C102-A$3*SIN(B102) 
E102: =C102-A$1*SIN(B102+A$2*C102) 
F102: =COS((D102+B102)/2) 
G102: =A$2*(E102+C102)-D102+B102 
H102: =-A$1*SIN((D102+B102)/2)-E102+C102
I102: =1/(1-A$4*F102) 
J102: =I102*(G102-A$2*H102) 
K102: =I102*(H102-A$2*F102*G102)

[6] 微分方程式の離散近似 （共著者の提案）
  http://pulsar.blog.ocn.ne.jp/IEICE/JHES2011-F4.pdf
[7] 広義の中点法 MP1, MP2, MP3 の比較
  http://pulsar.blog.ocn.ne.jp/IEICE/NOLTA2011-8037.pdf