Scilabで単振り子 その1 解析解との比較

これまで微分方程式による物理現象のモデル化(PDF)から運動学の章で紹介されている微分方程式をScilabを用いて計算してきました。
今回から4回は単振り子の微分方程式をScilabで計算します。

001_20130720182725.png

Fig.1: 単振り子の問題設定



単振り子の微分方程式と解析解


導出は次回に回しますが、単振り子の運動を表す微分方程式は、以下のように書くことが出来ます。

\frac{\mathrm{d}\theta}{\mathrm{d}t} = q
\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = - \frac{g}{L}\sin\theta

(θ: 振れ角, t: 時間, q: 角速度, g: 重力加速度, L: 糸の長さ)

この微分方程式は、解析的に解を求めることが出来ます。
例えば、初角度θ0から初角速度q0=0で振動を開始したとすると下記のようになります。

t(\theta) = - \sqrt{\frac{L}{2g}}\int^{\theta}_{\theta_0}\frac{\mathrm{d}\theta}{\sqrt{\cos\theta - \cos\theta_0}}

しかしながら、この積分は初等関数で表せる形に出来ないので、結局は数値積分をせざるを得ません。

今回は、これまでと同様にScilabの常微分方程式ソルバodeを用いて振りこの運動をシミュレーションするとともに、解析解に対する数値積分も行い、これら二つを比較します。

プログラミング


常微分方程式ソルバはode関数を使います。(参考:常微分方程式のタグが付いたエントリ)

解析解のほうの数値積分はintegrate関数を用いました。(参考:Scilabで数値積分: 固体の比熱,integrate - 求積法による積分)

// *** 解析解とソルバ解の共通部分 ***
// 初期振幅の入力
g = 9.8; // 重力加速度
l = g / (2 * %pi) ^ 2; // 糸の長さ
// 初期条件
th0 = input("th0 = "); // 初角度
q0 = 0; // 初角速度
X0 = [th0; q0];

// *** 常微分方程式ソルバによる解 ***
// 解くべき微分方程式の定義
function dx = pend(t,x)
dx(1) = x(2);
dx(2) = - g / l * sin(x(1));
endfunction
T = linspace(0,2,200); // 時間ベクトル
TH = ode(X0, 0, T, pend); // 常微分方程式ソルバ
plot(T,TH(1,:),'-r'); // プロット

// *** 解析解 ***
THanaly = linspace(- th0, th0, 20);
Tanaly = - sqrt(l / (2 * g)) * integrate('1 ./ sqrt(cos(th) - cos(th0))','th',th0,THanaly);
plot(Tanaly,THanaly,'og');

// *** グラフの体裁 ***
legend("o.d.e","integration",1);
xlabel("t[s]");
ylabel("$\theta \mathrm{[rad]}$");


計算結果


以下に、初角速度q0=0で初角度θ0を0.1ラジアン、2.9ラジアンの2通りの値で計算した結果を示します。

002_20130720182725.png

003_20130720182724.png

Fig.2-3: 厳密な単振り子の解。sinθ≒θの近似が成り立つときの周期が1秒となるように糸の長さLを調整している。振幅(初角度)が小さいとき(θ0=0.1; グラフ上)は周期が1秒となっているが、振幅(初角度)が大きいとき(θ0=2.9; グラフ下)は周期が1秒を大きく超えている。


関連エントリ




参考URL




付録


このエントリで使用したLTspiceのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。


tag: Scilab 常微分方程式 ode 数値積分 単振り子 

comment

Secret

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj常微分方程式モンテカルロ解析状態密度トランジスタodeDOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632AR6452A可変抵抗トランジスタ技術温度解析ブレッドボードI2C反強磁性確率論数値積分セミナーバンドギャップバンド構造偏微分方程式非線形方程式ソルバ熱設計絶縁ISO-I2Cカオス三端子レギュレータLM358GW近似マフィンティン半径A/DコンバータフォトカプラシュミットトリガLEDPC817C発振回路数値微分直流動作点解析サーボカレントミラーTL431アナログスイッチUSB74HC4053bzqltyVESTA補間電子負荷アセンブライジング模型BSch量子力学単振り子2ちゃんねるチョッパアンプLDA開発環境基本並進ベクトルFFT標準ロジックブラべ格子パラメトリック解析抵抗SMPMaxima失敗談ラプラス方程式繰り返し位相図スイッチト・キャパシタ熱伝導状態方程式キュリー温度gfortranコバルトTLP621不規則合金Quantum_ESPRESSO六方最密充填構造ランダムウォーク相対論ewidthスピン軌道相互作用FETQSGWVCAcygwinスレーターポーリング曲線GGA仮想結晶近似PWscfシュレディンガー方程式LM555ハーフメタル固有値問題NE555最小値ガイガー管QNAPUPS自動計測ダイヤモンドマントルTLP552格子比熱最適化MCU井戸型ポテンシャル最大値xcrysdenCIF条件分岐詰め回路フェルミ面差し込みグラフスーパーセルfsolveブラウン運動awk過渡解析起電力三角波第一原理計算FXA-7020ZRWriter509Ubuntuテスタ熱力学データロガーTLP521OpenMPubuntu平均場近似MAS830LトランスCK1026PIC16F785PGA2SC1815EAGLEノコギリ波負帰還安定性ナイキスト線図MBEOPA2277P-10フィルタCapSenseAACircuitLMC662文字列固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型電荷密度重積分SIC二相共存磁気モーメント不純物問題PWgui擬ポテンシャルゼーベック係数ZnOウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデルcif2cell入出力陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ