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)微分は1次近似の極限であり,極限をとる前でも基本的な特徴(例えば解軌道の曲がり方)が現れていることが多い.(2)Δt を変えても近似解があまり変わらなければほぼ収束しているとして扱えば気楽.
補足1: 微分方程式が線形の場合,ラプラス変換のs平面の原点近傍が中点法で離散近似したz変換の z =1 の近傍に比較的精度よく対応します.インパルス応答の標本値を保存する変換(標準z変換)があるにも拘わらず,そうでない変換(双1次変換)をフィルタの設計に用いるのはこのためです --- 虚軸上に極を持つような回路はフィルタに使えません.減衰振動の仕方も好ましく近似されていることが重要です.
補足2: 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
※コメント投稿者のブログIDはブログ作成者のみに通知されます