amazon：Excelで操る! ここまでできる科学技術計算 第2版 楽天： Excelで操る! ここまでできる科学技術計算 [ 神足 史人 ] |
 図１ 単振り子の説明 |
単振り子の運動方程式
単振り子の運動方程式は d2θ/dt2 = -(g/l)sinθ ・・・ (1) ｌ ：振り子の糸の長さ θ：糸が鉛直方向となす振れの角 ｔ ：経過時間 ｇ ：重力の加速度(9.8) と書ける。（振り子の質量(m)は振り子運動に影響しない） この方程式は非線形であって、解析解は楕円関数で表現される。 振れの角θが小さいとき、sinθ≒θとなり、(1)式は単振動の式 d2θ/dt2 = -(g/l)θ ・・・ (2) となる。 (2)式の解は θ＝2k・sin{√(g/l)・t} となり。振幅(k)、周期(T = ２π√(g/l))の単振動を表す。 (1)式の解析解は難しい楕円関数で表現されるが、ルンゲクッタ法に よる数値積分は容易である。 |
たった1日で即戦力になるExcelの教科書  入門者のExcel VBA―初めての人にベストな学び方 (ブルーバックス)  Excel VBAのプログラミングのツボとコツがゼッタイにわかる本 数値積分法の基礎と応用 参考資料：日本機械学会編 独修微分積分学改訂増補 数値解析  偏微分方程式の数値解法入門  数値積分法の基礎と応用  |
単振り子の１階連立微分方程式
単振り子の運動方程式（①式）は２階微分方程式である。 ルンゲクッタ法は１階連立微分方程式の解法である。そのため①式を１階連立微分方程式に変換する。 角速度(ω)は時間(sec)あたりの角度(ラジアン)変化である。よって、 dθ/dt = ω 単位換算 : 180度＝πラジアン です。 となる。 角速度（ω)を使うと、①式は dω/dt = -(g/l)sinθ dθ/dt = ω の１階連立微分方程式で表現できる。また関数F、Gを使い、 dω/dt = F( t , ω, θ) F( t , ω, θ) = -(g/l)sinθ dθ/dt = G( t , ω, θ) G( t , ω, θ)= ω と表現できる。 |
４次のルンゲクッタ法
次の１階連立微分方程式で説明する。 dω/dt = F( t , ω, θ) dθ/dt = G( t , ω, θ) F,Gは時間tと変数 ω,θの関数 初期値 t=t0 の時、ω=ω0,θ=θ0 とすると、Δt後の変数値は下式となる。 k1 = Δt・F( t0, ω0, θ0) m1 = Δt・G( t0, ω0, θ0) k2 = Δt・F( t0+Δt/2,ω0+k1/2, θ0+m1/2) m2 = Δt・G( t0+Δt/2, ω0+k1/2, θ0+m1/2) k3 = Δt・F( t0+Δt/2, ω0+k2/2, θ0+m2/2) m3 = Δt・G( t0+Δt/2, ω0+k2/2, θ0+m2/2) k4 = Δt・F( t0+Δt, ω0+k3, θ0+m3) m4 = Δt・G( t0+Δt, ω0+k3, θ0+m3) k = (k1+2・K2+2・k3+k4)/6 m = (m1+2・m2+2・m3+m4)/6 Δt後の変数値ω1,φ1は ω1 = ω0 + k θ1 = θ0 + m となる。 この手続きを繰り返せば、離散的だが、すべての時間での変数値が求まる。 ちなみに、Δtは0.01秒のような微小値を設定する。 Δtが大きすぎると解が発散する場合がある。 計算量が多いので、手計算は困難です。EXCELなどを使って計算しましょう。 |
エクセルを用いたルンゲクッタ法
ルンゲクッタ法を用い単振り子の2階連立微分方程式を解いて見ましょう。 振り子の糸の長さ(l) ＝ 0.5ｍ とします。 振れの角度(θ)＝ 1ラジアン = 180/π度 = 57.3度からスタートします。 すなわち、初期値は ｔ0＝０秒 ω0=0ラジアン/秒 θ0=1ラジアン となります。 エクセルを用いたルンゲクッタ法 を参照し、エクセルシートに計算式を作り、繰り返し計算のためのVBAを作って 下さい。  図２ 単振り子の位置、角速度を解析するためのエクセルシート 計算の繰り返し処理のため、ＶＢＡです。  図３ ルンゲクッタ法のための計算繰り返し用VBA ルンゲクッタ法による微分方程式の計算結果を図４と図５のグラフにしました。  図４  図５ |
振り子振動の周期と等時性
単振動の周期の理論解は(3)式となる。 T0 = ２π√(ｌ/g) ・・・・・ (3) 周期とは、振り子が１回振動するのに要する時間(秒)である。 振り子の糸の長さ(l) ＝ 0.5ｍ の場合、(3)式より T0 = 1.4192秒 となる。 振れの角度(θ)が大きい場合、周期は近似的に(4)式となり、そのグラフ は図6となる。 T = T0(1+θ2/16) ・・・・・ (4)  図６ 振れの角度(θ)=1の場合、 T = 1.5079秒 となり、ルンゲクッタ法による解析結果と良く一致している。 ふれの角度が１ラジアン=57°以下ならば振り子の周期はほとんど振幅に 無関係である。 これは振り子の等時性と呼ばれる性質である。 |
振り子に大きな初速（角速度）を設定すると？
単振り子に大きな初速（角速度)を設定し、運動方程式を解くとどうなるでしょう。 初期角速度(ω)＝ 9ラジアン/秒 を設定してみましょう。 すなわち、初期値は ｔ0＝０秒 ω0=9ラジアン/秒 θ0=0ラジアン とする。  図7 大きな初速度を振り子に設定した場合の解析結果とそのグラフ 初速（角速度）が大きい場合、単振り子は振動せず、回転しています。（振り子ではありません） |
参考文献
数値積分法の基礎と応用
日本機械学会編
Runge-Kutta-Fehlberg法
新潟大学工学部化学システム工学科からの依頼により、エクセルを用いたRunge-Kutta-Fehlberg法を製作しました。Runge-Kutta-Fehlberg法は誤差をチェックし、時間刻み幅を変更しながら計算するため、通常のルンゲクッタ法と比較して非常に精度が良い方法です。通常のルンゲクッタ法では発散してしまう温度境界層問題も解けます。Runge-Kutta-Fehlberg法によるVBAは化学工学(反応速度式などに適応)のみならず機械工学にもたいへん有用と思われます。この常微分方程式解法シート(ITOシート)は新潟大学の化学工学資料のページからダウンロードできます。 参考 ： 化学工学資料のページ Excelで気軽に化学工学 
エクセルを用いたルンゲクッタ法による多質点系バネマスモデルの振動解析 上記振動解析において、４次のルンゲクッタ法とRunge-Kutta-Fehlberg法を比較しました。Runge-Kutta-Fehlberg法はエネルギー保存則(運動エネルギー＋バネのエネルギー)を満たしており、ひじょうに精度の良い解法であることがわかる。また４次のルンゲクッタ法より精度が良いにもかかわらずその演算時間が格段に小さく、性能が優れていることが解る。 エクセルを用いたルンゲクッタフェールベルグ法によるクロソイド曲線(緩和曲線) 計算が難しいといわれているクロソイド曲線も容易に算出できる。 ルンゲクッタフェールベルグ法により電気回路方程式を解く 電気回路にも適用できる。 Excelを用いてラグランジェポイントにおける物体の運動方程式を解く (決して宇宙人には教えないで下さい) ラグランジェポイントでは物体は奇妙な運動を行う。 |
エクセルを用いたルンゲクッタ法による多質点系モデルの振動解析
ルンゲクッタフェールベルグ法により電気回路方程式を解く
Excelを用いてラグランジュポイントにおける物体の運動方程式を解く
(決して宇宙人には教えないで下さい)
エクセルを用いたルンゲクッタフェールベルグ法によるクロソイド曲線(緩和曲線)
エクセルのソルバーを用いた最小作用の原理に基づく質点の軌道計算
大発見、エクセルのソルバーで質点の軌道計算ができる。
質点の軌道に関する理論解と適当に設定した軌道
Excelを用いた地球温暖化計算
EXCELを用いた科学技術計算
EXCELのソルバーを用いた科学技術計算
EXCELを用いた構造力学(材料力学)計算
〒６７３－００３６ 兵庫県明石市松江６２－１４
(有)ゴッドフット企画
感想をお聞かせ下さい
ＪＲ福知山線 脱線事故シミュレーション
ルンゲクッタ法で脱線事故シミュレーション
（Yahoo Japan 掲載）
Copyright (C) GODFOOT 2007 All Rights Reserved.