VCSSL > 日本語トップ > コードアーカイブ > シミュレーション

[ 前へ | 目次 | 次へ ]

Now Loading...

PC （※スマートフォンでは動きません）でダウンロードし、ZIPファイルを右クリックメニューから展開して、できたフォルダ内の「 VCSSL.bat（バッチファイル）」をダブルクリックすると起動します。 Linux等では「 VCSSL.jar 」をコマンド実行してください。
» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら

二重振り子のシミュレーション

解析力学のラグランジュ方程式をオイラー法で数値的に解く事で、二重振り子の挙動をアニメーション表示するVCSSLプログラムです。

- 目次 -

使用方法
題材解説
コード解説
ライセンス

使用方法

ダウンロードと展開（解凍）

まず、PC（スマホは未対応）で上の画面の「ダウンロード」ボタンを押してください。するとZIP形式で圧縮されたファイルがダウンロードされます。

Windows をご使用の方は、ここでまずZIPファイルを右クリックし、「プロパティ」を選んで開かれる画面で、下の方にあるセキュリティ項目の「許可する」にチェックを入れて「OK」で閉じてください。これを行わないと、ZIP展開やソフト起動時に、警告メッセージが出て展開完了/起動できない場合があります。

その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開（解凍）してください。展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…

プログラムの起動

Windows をご使用の場合

上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください：

VCSSL__ダブルクリックでプログラム実行.bat

もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。

正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。

» うまく起動できずにエラーになってしまう場合は…

Linux 等をご使用の場合

ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください：

java -jar VCSSL.jar
（プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行）

» javaコマンドが使用できない等のエラーが表示される場合は…

操作方法

このプログラムは、アルゴリズム紹介のための簡易プログラムであるため、操作機能などは付いていません。パラメータ等の設定は、コード中の記述を直接変更してください。

題材解説

二重振り子のモデル

このシミュレーションのモデルは、普通の振り子の先に、もう一つの振り子がぶら下がった、二重振り子と呼ばれるものです。

二重振り子は、中心の支点から1つめの振り子がぶら下がり、そこからさらに2つめの振り子がぶら下がっています。

二重振り子の図

青いボールが1つめの振り子で、その鉛直線からの角度がθ1 です。そこからぶら下がる赤いボールが2つめの振り子で、その鉛直線からの角度がθ2 です。

図のように、1つめの振り子が鉛直線となす角度をθ1、2つめの振り子が円直線となす角度をθ2とします。

なお、θの時間微分、つまり回転のスピードは、上にドットを付けてθ・と表記します。さらにその時間微分、つまり回転の加速度は、ドットを2つ付けてθ・・と表記します。

また、2つの振り子のバーの長さをそれぞれL1とL2、ボールの重さをそれぞれM1とM2と表記します。

ニュートン方程式で素直に扱うと難しいので、解析力学で扱う

このシミュレーションでは、一般によく使用されるような、ニュートンの運動方程式は使用しません。二重振り子のような複雑なモデルをニュートン方程式で扱うと、かなり難しくなってしまうためです。

そこでこのシミュレーションでは、ニュートン方程式の代わりに、「解析力学」という分野の、少し高度な方法を使用します。解析力学の手法を用いれば、複雑なモデルを、ニュートン方程式よりもシンプルに扱う事ができます。

解析力学を用いた場合、モデルのエネルギーを表す数式さえ分かれば、そのモデルの自由度変数（「動き」を指定する変数、二重振り子なら角度）に関する加速度を、直接的に求める事ができます。

二重振り子の場合、2つの振り子の角度に関する加速度さえ分かれば、あとはオイラー法などですぐにシミュレーションできますね。振り子以外でも同様です。自由度変数に関する加速度さえ分かれば、そのモデルの動きは解析できます。自由度変数以外は、そもそも動かないのですから。

ラグランジアンとは

まず求めるべきなのが、ラグランジアンと呼ばれる特殊な物理量です。

ラグランジアンそのものに関する詳しい説明は解析力学の教科書に譲るとして、要するに、これは「運動エネルギー - 位置エネルギー」で求まる量です。

$L=T-U$

上のLがラグランジアン、Tが運動エネルギー、Uが位置エネルギーです。なおUに関しては、もう少し厳密には「ポテンシャルエネルギー」なのですが、日常的なモデルでは位置エネルギーと言っても良いでしょう。

運動エネルギー T を求める

まずは運動エネルギー T を求めましょう。

各ボールの速度V1とV2は、振り子の運動を垂直/水平に分けて考えれば：

$V_{1}=(L_{1}\dot\theta_{1}\cos\theta_{1})^{2}+(L_{1}\dot\theta_{1}\sin\theta_{1})^{2}$
$=L_{1}^2\dot\theta_{1}^2$

$V_{2}=(L_{1}\dot\theta_{1}\cos\theta_{1}+L_{2}\dot\theta_{2}\cos\theta_{2})^{2}+(L_{1}\dot\theta_{1}\sin\theta_{1}+L_{2}\dot\theta_{2}\sin\theta_{2})^{2}$
$=L_{1}^2\dot\theta_{1}^2+L_{2}^2\dot\theta_{2}^2+2L_{1}L_{2}\dot\theta_{1}\dot\theta_{2}\cos(\theta_1-\theta_2)$

従って、このモデルの運動エネルギー T は

$T=\frac{1}{2}M_{1}V_{1}^{2}+\frac{1}{2}M_{2}V_{2}^{2}$
$=\frac{1}{2}(M_{1}+M_{2})L_{1}^{2}\dot\theta_{1}^{2}+\frac{1}{2}M_{2}L_{2}^{2}\dot\theta_{2}^{2}+M_{2}L_{1}L_{2}\dot\theta_{1}\dot\theta_{2}\cos(\theta_1-\theta_2)$

と求まりました。

位置エネルギー U を求める

続いて、位置エネルギー U を求めましょう。

位置エネルギーの基準となる高さは、どこにとっても良いので、ここでは1つめの振り子の支点の高さにとりましょう。すると、ボールが支点よりも下にある時に、位置エネルギーがマイナスになりますが、特に問題はありません。

この場合、このモデルの位置エネルギー U は、支点から鉛直上向きにy軸をとって：

$U=M_{1}gy_{1}+M_{2}gy_{2}$
$=-M_{1}gL_{1}\cos\theta_{1}-M_{2}g(L_{1}\cos\theta_{1}+L_{2}\cos\theta_2)$

と求まりました。ここでgは重力加速度です。

ラグランジアンとラグランジュ方程式

運動エネルギーと位置エネルギーが求まったので、それらの差として、ラグランジアンの数式を以下のように求める事ができます。

$L=\frac{1}{2}(M_{1}+M_{2})L_{1}^{2}\dot\theta_{1}^{2}+\frac{1}{2}M_{2}L_{2}^{2}\dot\theta_{2}^{2}+M_{2}L_{1}L_{2}\dot\theta_{1}\dot\theta_{2}\cos(\theta_1-\theta_2)$
$-M_{1}gL_{1}\cos\theta_{1}-M_{2}g(L_{1}\cos\theta_{1}+L_{2}\cos\theta_2)$

さて、ただラグランジアンが求まっただけでは、まだ運動の様子は分かりません。「力から運動の様子を求める」ニュートンの運動方程式のように、「ラグランジアンから運動の様子を求める」ための基礎方程式が必要です。この基礎方程式となるのが、オイラー = ラグランジュ方程式、通称ラグランジュ方程式です。

このモデルの場合、ラグランジュ方程式は以下のように2つ与えられます：

$\frac{d}{dt}\frac{\partial L}{\partial\dot\theta_{1}}-\frac{\partial L}{\partial\theta_{1}}=0$

$\frac{d}{dt}\frac{\partial L}{\partial\dot\theta_{2}}-\frac{\partial L}{\partial\theta_{2}}=0$

分野によっては、∂L/∂θといった表記は見慣れないものかもしれません。これは偏微分と呼ばれるもので、θ以外の変数をすべて定数と見なす、やや強引な微分を意味します。θ・やθ・・も定数と見なします。つまり、特に何も考えず、θだけに着目して微分すれば良いのです。

そして∂L/∂θ・では、今度はθ・のみを変数と見なし、θやθ・・などはすべて定数と見なして微分します。

そうしてラグランジュ方程式中の微分を全て計算すると、以下の2つの等式が得られます：

$(M_{1}+M_{2})L_{1}^{2}\ddot\theta_{1}+M_{2}L_{1}L_{2}\cos(\theta_{1}-\theta_{2})\ddot\theta_{2}$
$=-M_{2}L_1L_2\dot\theta_{2}^{2}\sin(\theta_{1}-\theta_{2})-(M_{1}+M_{2})gL_{1}\sin\theta_{1}$ 　　(1)

$M_{2}L_{1}L_{2}\cos(\theta_{1}-\theta_{2})\ddot\theta_{1}+M_{2}L_{2}^{2}\ddot\theta_{2}$
$=M_{2}L_1L_2\dot\theta_{1}^{2}\sin(\theta_{1}-\theta_{2})-M_{2}gL_{2}\sin\theta_{2}$ 　　(2)

こうして無事、ラグランジアンから、ラグランジュ方程式を用いて、二重振り子の運動を支配する方程式が導かれました。

行列化

さて、上で求めたラグランジュ方程式は、θ1・・とθ2・・の連立方程式になっているのが少々厄介です。実際にオイラー法などで運動をシミュレーションするには、この連立方程式を解き、以下のような形に分解する必要があります：

$\ddot\theta_{1}=...$ 　　(3)

$\ddot\theta_{2}=...$ 　　(4)

問題は、実際にどうやってこの連立方程式を解くかです。

素直にやるなら、(1)式と(2)式に何かしらの値をかけたり、互いに引いたり足したりして、まず θ1・・かθ2・・のどちらかを消去し、そして得られた値を他方に代入する事になります。

しかし、問題が拡張されて振り子の数が増え、例えば百重振り子や千重振り子などになってくると、そのようなやり方は通用しなくなります。1千階の連立方程式などを手で解くのは非現実的だからです。

そこで今回は応用性を考え、やや大げさですが、行列を使って解いてみましょう。このやり方なら、例え千重振り子にでも応用できるはずです。

まず、(1)式と(2)式の連立方程式を、以下のように行列形式に再記述します：

$($	$A_{1}$	$B$	$)$
	$B$	$A_{2}$

$($	$\ddot\theta_{1}$	$)$
	$\ddot\theta_{2}$

$=$

$($	$D1$	$)$
	$D2$

　(5)

ここで行列中の各文字は、以下のように置きました。

$A_{1}=(M_{1}+M_{2})L_{1}^{2}$ 　(6)

$A_{2}=M_{2}L_{2}^{2}$ 　(7)

$B=M_{2}L_{1}L_{2}\cos(\theta_{1}-\theta_{2})$ 　(8)

$D_{1}=-M_{2}L_{1}L_{2}\dot\theta_{2}\sin(\theta_{1}-\theta_{2})-(M_{1}+M_{2})gL_{1}\sin(\theta_{1})$ 　(9)

$D_{1}=M_{2}L_{1}L_{2}\dot\theta_{1}\sin(\theta_{1}-\theta_{2})-M_{2}gL_{2}\sin(\theta_{2})$ 　(10)

逆行列を用いて、行列式を解く

さて、(5)のように行列式にしてしまえば、あとは簡単です。左辺にある行列の逆行列を求めて、両辺に左からかけるだけです。実際にやってみましょう。

2次元行列の逆行列を求めるのは、高校数学でも習う、以下の公式を使えば一瞬です：

$($	$a$	$b$	$)$	$-1$
	$c$	$d$

$=$

$\frac{1}{ad-bc}$	$($	$d$	$-b$	$)$
		$-a$	$c$

従って(5)式左辺にある行列の逆行列は：

$($	$A_{1}$	$B$	$)$	$-1$
	$B$	$A_{2}$

$=$

$\frac{1}{A_{1}A_{2}-B^{2}}$	$($	$A_{2}}$	$-B$	$)$
		$-B$	$A_{1}$

と求まります。

実は2次元だけでなく、3次元、4次元、さらに大きな次元の行列についても、逆行列を求める公式が存在します。さらに、コンピューターで数値的に逆行列を求めるアルゴリズムも完成されており、行列演算ライブラリに当然のように入っています。大きな次元の場合は、それらを利用しましょう。

さて、この逆行列を(5)式の両辺に左からかけて：

$($	$A_{1}$	$B$	$)$	$-1$
	$B$	$A_{2}$

$($	$A_{1}$	$B$	$)$
	$B$	$A_{2}$

$($	$\ddot\theta_{1}$	$)$
	$\ddot\theta_{2}$

$=$

$\frac{1}{A_{1}A_{2}-B^{2}}$	$($	$A_{2}}$	$-B$	$)$
		$-B$	$A_{1}$

$($	$D1$	$)$
	$D2$

各成分を実際に計算すると：

$($	$1$	$0$	$)$
	$0$	$1$

$($	$\ddot\theta_{1}$	$)$
	$\ddot\theta_{2}$

$=$

$\frac{1}{A_{1}A_{2}-B^{2}}$	$($	$A_{2}D_{1}-BD_{2}$	$)$
		$A_{1}D_{2}-BD_{1}$

よって解は：

$\ddot\theta_{1}=\frac{A_{2}D_{1}-BD_{2}}{A_{1}A_{2}-B^{2}}$ 　(11)

$\ddot\theta_{2}=\frac{A_{1}D_{2}-BD_{1}}{A_{1}A_{2}-B^{2}}$ 　(12)

と、このように求まりました。なおA、B、Dの値に関しては(6)～(10)式を参照してください。

オイラー法

最後に、今回のシミュレーションに用いる技法である、オイラー法についてです。ここではあまり詳しく触れませんが、イメージだけでも確認しておきましょう。

まず、微分を以下のように近似します。

$\ddot\theta=\frac{d\dot\theta}{dt}\approx\frac{\Delta\dot\theta}{\Delta t}$ 　(13)

$\dot\theta=\frac{d\theta}{dt}\approx\frac{\Delta\theta}{\Delta t}$ 　(14)

ここで Δt は「時間刻み」や「タイムステップ」と呼ばれる量であり、いってみればアニメーションのコマ撮り間隔のようなものです。この値を細かく設定するほど、シミュレーションの精度が向上する反面、処理速度が遅くなります。

これらの近似式の両辺にΔtをかけると、Δtの微小時間におけるθ・とθの変位が、以下のように求まります：

$\Delta\dot\theta=\ddot\theta\Delta t$

$\Delta\theta=\dot\theta\Delta t$

よって、現在時刻の状態から、次の時刻の状態を求める漸化式が、以下のように求まります。

$\dot\theta_{1}^{{\rm next}}=\dot\theta_{1}+\Delta\dot\theta_{1}=\dot\theta_{1}+\ddot\theta_{1}\Delta t=\dot\theta_{1}+\frac{A_{2}D_{1}-BD_{2}}{A_{1}A_{2}-B^2}\Delta t$ 　(15)

$\dot\theta_{2}^{{\rm next}}=\dot\theta_{2}+\Delta\dot\theta_{2}=\dot\theta+\ddot\theta_{2}\Delta t=\dot\theta_{2}+\frac{A_{1}D_{2}-BD_{1}}{A_{1}A_{2}-B^2}\Delta t$ 　(16)

$\theta_{1}^{{\rm next}}=\theta+\Delta\theta_{1}=\theta_{1}+\dot\theta_{1}\Delta t$ 　(17)

$\theta_{2}^{{\rm next}}=\theta_{2}+\Delta\theta_{2}=\theta+\dot\theta_{2}\Delta t$ 　(18)

これで、シミュレーションに必要な計算は全て終了しました。あとは、時刻をΔt進めるごとに、(15)～(18)式をひたすら繰り返し処理し続ければ良いだけです。

なお、BとDの値もθやθ・に依存しているため、これらの値も毎時刻、(6)～(10)式を用いて再計算してやる必要があります。

コード解説

コーディング

ここまでで、必要な計算とアルゴリズムの解説は終わったので、ここからはそれをコーディングし、シミュレーションプログラムとして仕上げましょう。

最初に、少し長いですが、コード全体を掲載します。

DuplexPendulum.vcssl

以上の通りです。それでは、細部を見ていきましょう。

メインループ型の処理フロー

今回のシミュレーションでは、メインループ型の処理フローを採用します。これは、プログラムの中核部分が一種の無限ループ（メインループ）になっており、プログラムの開始から終了まで、ずっとループが回り続けているようなフローです。

典型的なメインループ型処理フローは、下図のようになります。

典型的な、メインループ型処理フローの例

まずmain関数の最初の部分でinit関数がコールされ、初期化などの処理が行われます。そしてメインループに入り、プログラム終了まで、描画（paint関数）、更新（update関数）、待機（sleep関数）を繰り返し続けます。ウィンドウが閉じられ、onWindowClose関数が呼ばれると、メインループを脱出し、main関数の最後でexit関数をコールして、プログラムを終了します。

上図のように、メインループ内で毎回、画面描画と情報更新を繰り返し続けます。

通常、メインループ型処理フローでは、ユーザーからの操作があった場合に、それに対応する処理もメインループ内で行います。

しかし今回はあまりインタラクティブなものは作らないので、ウィンドウを閉じる操作があった場合のみ、プログラムを終了させる処理を行います。

テンプレート

さて、メインループ型処理フローを採用する場合、コードの大まかな基本構造は、ある程度決まりきったものになります。それを毎回ゼロからコーディングするのは面倒なので、あらかじめメインループ型処理フローの基本構造が実装されたテンプレートを使用しましょう。

- メインループ型テンプレート -

// ヘッダ領域
// ここでライブラリのインポートや、グローバル変数の定義など

bool loopState = true; // メインループを脱出させるための変数
int loopWait = 24; // メインループ1周ごとの待機時間（ミリ秒）

// プログラムはここから開始する
void main( string args[] ){
    init(); // 初期化
    // メインループ
    while( loopState ){
        paint(); // 描画
        update(); // 更新
        sleep( loopWait ); // 待機
    }
    exit(); // メインループを抜けたらプログラムを終了
}

// 起動後の初期化処理を行う関数
void init(){
    // ここに初期化処理を記述
}

// 毎周期の更新処理を行う関数
void update(){
    // ここに更新処理を記述
}

// 毎周期の画面描画を行う関数
void paint(){
    // ここに描画処理を記述
}

// イベントハンドラ、ウィンドウが閉じられた際にコールされる
void onWindowClose( int id ){
    loopState = false; // メインループを脱出、プログラム終了
}

このテンプレートのコードは、「何もしないメインループ型のプログラム」となっています。何もしないだけで、ちゃんと実行するとメインループが回り続け、ウィンドウを閉じると終了します。

テンプレートの各部に、必要に応じて目的の処理を追記していく事によって、シミュレーションプログラムへと仕上げていきます。

ヘッダ領域の記述

まずは、プログラムの先頭、ヘッダ領域を記述しましょう。ここでは、文字コードの指定やライブラリのインポート、グローバル変数の宣言などを行います。

最初の行には、coding宣言を用いて文字コードを指定します。これは必須ではありませんが、指定しておくと文字化けによる不具合を防ぐ事ができます。使用可能な文字コードは「Shift_JIS」または「UTF-8」です。

続いてimport宣言で、必要なライブラリをインポートします。今回は、GUIやグラフィックス関連、および数学関数関連の、計4つの標準ライブラリを使用します。

coding Shift_JIS;

import GUI;
import Graphics;
import Graphics2D;
import Math;

グローバル変数の宣言

import宣言よりも下では、必要なグローバル変数を宣言します。今回は以下の変数をグローバル変数として宣言しました。

// GUI用の変数
int WIDTH = 600; // ウィンドウの幅
int HEIGHT = 600; // ウィンドウの高さ
int window; // ウィンドウのID
int displayLabel; // 表示ラベルのID

// グラフィックス用の変数
int PIXEL_PER_LENGTH = 100; // 長さとピクセルの変換係数
int graphics; // グラフィックスのID
int renderer; // レンダラーのID

// モデルの性質を表すパラメータ変数
float gravityAccel = 9.8; // 重力加速度
float deltaT = 0.00001; // 時間刻み
float length1 = 1.4; // 振り子1の長さ
float length2 = 1.2; // 振り子2の長さ
float mass1 = 1.2; // 振り子1の重さ
float mass2 = 0.1; // 振り子2の重さ

// モデルの状態を表す自由度変数
float theta1 = 3.1; // 振り子1の角度
float theta2 = -2.0; // 振り子2の角度
float theta1Dot = 0.0; // 振り子1の角速度
float theta2Dot = 0.0; // 振り子2の角速度

// メインループ制御用の変数
bool loopState = true; // メインループを脱出させるための変数
int loopWait = 24; // メインループ1周ごとの待機時間（ミリ秒）

「GUI用の変数」の箇所では、ウィンドウなどのGUIコンポーネントを生成した際の、IDを格納する変数などを指定します。また、ウィンドウのサイズなども宣言しておきます。

「グラフィックス用の変数」の箇所では、グラフィックスデータや2DCGレンダラー（描画エンジン）を生成した際の、IDを格納する変数などを指定します。

「モデルの性質を表すパラメータ変数」の箇所では、二重振り子の長さや重さなど、静的なパラメータを宣言し、値を定義しています。

「モデルの状態を表す自由度変数」の箇所では、二重振り子の動的な量、つまり角度や角速度などを定義しています。この初期値を変える事よって、様々な振り子の運動が楽しめるでしょう。

「メインループ制御用の変数」の箇所は、テンプレートに最初からありますね。これらの変数は特に触る必要はありません。

中身が空の関数は、目的に応じた処理を独自に実装する

さて、テンプレートには、中身が空の関数がいくつか宣言されています。この中に、今回の二重振り子など、それぞれの用途に応じた処理を実装していきます。それぞれの関数の概要は下記の通りです。

init関数 … 最初に一度だけ実行される。初期化などを行う。
paint関数 … メインループが1周するごとに毎回呼ばれる。画面の描画処理などを行う。
update関数 … メインループが1周するごとに毎回呼ばれる。情報の更新処理などを行う。

これら3つの関数に、今回のシミュレーションの処理を実装していきます。

init関数 ― 初期化処理の実装

まずは、初期化処理です。これはinit関数の中に実装します。「初期化」というのは、簡単に言えば、準備のようなものです。

今回は、ここでグラフィックスデータと2DCGレンダラー（描画エンジン）の生成と、表示画面の生成を行います。

// 起動後の初期化処理を行う関数
void init(){
    // グラフィックスデータを生成
    graphics = newGraphics();
    // 上で生成したgraphicsに描き込む2DCGレンダラーを生成
    renderer = newGraphics2DRenderer( WIDTH, HEIGHT, graphics );

    // ウィンドウを生成
    window = newWindow( 0, 0, WIDTH, HEIGHT+20, "Duplex Pendulum" );
    // graphicsの内容を表示するグラフィックスラベルを生成
    displayLabel = newImageLabel( 0, 0, WIDTH, HEIGHT, graphics );
    // ラベルをウィンドウに配置
    mountComponent( displayLabel, window );
}

これで初期化処理の実装は終了です。

この時点でプログラムを実行すると、空白のウィンドウが表示されます。まだ画面描画の処理を実装していないからです。

paint関数 ― 描画処理の実装

それでは、引き続き描画処理を実装しましょう。これはpaint関数の中で行います。

// 毎周期の画面描画を行う関数
void paint(){
    // 背景を白で塗りつぶし
    setDrawColor( renderer, 255, 255, 255, 255 );
    drawRect( renderer, 0, 0, WIDTH, HEIGHT, true );

    // 振り子のボールの座標を求める
    int centerX = WIDTH / 2;
    int centerY = HEIGHT / 2;
    int point1X = PIXEL_PER_LENGTH * length1 * sin( theta1 ) + centerX;
    int point1Y = PIXEL_PER_LENGTH * length1 * cos( theta1 ) + centerY;
    int point2X = PIXEL_PER_LENGTH * length2 * sin( theta2 ) + point1X;
    int point2Y = PIXEL_PER_LENGTH * length2 * cos( theta2 ) + point1Y;

    // 黒色で、振り子のバーを描画
    setDrawColor( renderer, 0, 0, 0, 255 );
    drawLine( renderer, centerX, centerY, point1X, point1Y );
    drawLine( renderer, point1X, point1Y, point2X, point2Y );

    // それぞれの色で、振り子のボールを描画
    int r = 10;
    setDrawColor( renderer, 0, 255, 0, 255 );
    drawOval( renderer, centerX-r, centerY-r, 2*r, 2*r, true );

    setDrawColor( renderer, 0, 0, 255, 255 );
    drawOval( renderer, point1X-r, point1Y-r, 2*r, 2*r, true );

    setDrawColor( renderer, 255, 0, 0, 255 );
    drawOval( renderer, point2X-r, point2Y-r, 2*r, 2*r, true );

    // GUIコンポーネントの再描画
    paintComponent( displayLabel );
    paintComponent( window );
}

これで描画処理の実装は終了です。この時点でプログラムを実行すると、画面に二重振り子の2DCGが描画されます。しかし、更新処理を実装していないため、まだ動きません。

update関数 ― 更新処理の実装

最後に、更新処理の実装です。

更新処理では、モデルの状態を表す変数の値を更新します。モデルの状態は、リアルタイムで画面に描画され続けていますから、これでようやくモデルが動く様子が見えるわけです。

二重振り子の運動については、すでに上で解析力学を用いた計算を行い、(15)～(18)式として求まりました。その結果をもう一度見ておきましょう。オイラー法で二重振り子の状態を更新する漸化式は：

$\dot\theta_{1}^{{\rm next}}=\dot\theta_{1}+\frac{A_{2}D_{1}-BD_{2}}{A_{1}A_{2}-B^2}\Delta t$ 　(15)

$\dot\theta_{2}^{{\rm next}}=\dot\theta_{2}+\frac{A_{1}D_{2}-BD_{1}}{A_{1}A_{2}-B^2}\Delta t$ 　(16)

$\theta_{1}^{{\rm next}}=\theta_{1}+\dot\theta_{1}\Delta t$ 　(17)

$\theta_{2}^{{\rm next}}=\theta+\dot\theta_{2}\Delta t$ 　(18)

ここでA、B、Dの値は：

$A_{1}=(M_{1}+M_{2})L_{1}^{2}$ 　(6)

$A_{2}=M_{2}L_{2}^{2}$ 　(7)

$B=M_{2}L_{1}L_{2}\cos(\theta_{1}-\theta_{2})$ 　(8)

$D_{1}=-M_{2}L_{1}L_{2}\dot\theta_{2}\sin(\theta_{1}-\theta_{2})-(M_{1}+M_{2})gL_{1}\sin(\theta_{1})$ 　(9)

$D_{1}=M_{2}L_{1}L_{2}\dot\theta_{1}\sin(\theta_{1}-\theta_{2})-M_{2}gL_{2}\sin(\theta_{2})$ 　(10)

でした。ここでは、これらの漸化式を、そのまま素直に実装するだけです。

// 毎周期の更新処理を行う関数
void update(){
    float theta1DotDot, theta2DotDot;
    float a1, a2, b, d1, d2;

    // パラメータA1,A2の計算
    a1 = (mass1 + mass2) * length1 * length1;
    a2 = mass2 * length2 * length2;

    // アニメーションの1フレームあたり、オイラー法を5000ほど回す
    for( int skip=0; skip<5000; skip++ ){
        // パラメータB,D1,D2の計算
        b = mass2 * length1 * length2 * cos( theta1 - theta2 );
        d1 = -mass2 * length1 * length2 * theta2Dot * theta2Dot * sin( theta1 - theta2 ) - ( mass1 + mass2 ) * gravityAccel * length1 * sin( theta1 );
        d2 = mass2 * length1 * length2 * theta1Dot * theta1Dot * sin( theta1 - theta2 ) - mass2 * gravityAccel * length2 * sin( theta2 );

        // オイラー法で、振り子をdeltaTだけ運動させる
        theta1DotDot = ( a2*d1 - b*d2 ) / ( a1*a2 - b*b );
        theta2DotDot = ( a1*d2 - b*d1 ) / ( a1*a2 - b*b );
        theta1Dot += theta1DotDot * deltaT;
        theta2Dot += theta2DotDot * deltaT;
        theta1 += theta1Dot * deltaT;
        theta2 += theta2Dot * deltaT;
    }
}

以上で、二重振り子のシミュレーションプログラムは完成です。

ライセンス

このVCSSL/Vnanoコード（拡張子が「.vcssl」や「.vnano」のファイル）は実質的な著作権フリー（パブリックドメイン) である CC0 の状態で公開しています※。記事中にC言語/C++/Java言語などでのサンプルコードが掲載されいてる場合は、それらについても同様です。そのままでのご利用はもちろん、改造や流用などもご自由に行ってください。

※ ただし、このコードの配布フォルダ内には、ダウンロード後すぐに実行できるように、 VCSSLの実行環境も同梱されており、そのライセンス文書は「 License 」フォルダ内に同梱されています（要約すると、商用・非商用問わず自由に使用できますが、使用の結果に対して開発元は一切の責任を負いません、といった具合の内容です）。配布フォルダ内の各構成物の一覧やライセンスについては「 ReadMe_使用方法_必ずお読みください.txt 」をご参照ください。

※ Vnano の実行環境については、別途スクリプトエンジンのソースコードも一般公開しており、何らかのソフトウェア内に組み込んでご利用いただく事も可能です。詳細はこちらをご参照ください。

この記事中の商標などについて

OracleとJavaは、Oracle Corporation 及びその子会社、関連会社の米国及びその他の国における登録商標です。文中の社名、商品名等は各社の商標または登録商標である場合があります。
Windows は、米国 Microsoft Corporation の米国およびその他の国における登録商標です。この記事は独立著作物であり、Microsoft Corporation と関連のある、もしくはスポンサーを受けるものではありません。
Linux は、Linus Torvalds 氏の米国およびその他の国における商標または登録商標です。
その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。

[ 前へ | 目次 | 次へ ]

VCSSL > 日本語トップ > コードアーカイブ > シミュレーション

Vnano版 \| ローレンツ方程式を数値的に解くプログラム
	ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。
波の干渉（面上の円形波）のアニメーション表示
	面上の円形波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。
円形波のアニメーション表示
	振幅・波長・周期をスライダ―で操作しながら、円形波のグラフをアニメーションで見られるプログラムです。
波の干渉（線上の正弦波）のアニメーション表示
	線上（1次元の）の正弦波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。
正弦波のアニメーション表示
	振幅・波長・周期をスライダ―で操作しながら、正弦波のグラフをアニメーションで見られるプログラムです。
凹レンズを通過する波のシミュレーション
	凹レンズ形状の高密度媒質を通過する、波のシミュレーションです。
凸レンズを通過する波のシミュレーション
	凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。
乱雑な密度分布における波のシミュレーション
	密度分布が乱雑な媒質中における、波の伝播のシミュレーションです。
ローレンツアトラクタ（ファイル出力版）
	4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。
波の屈折のシミュレーション
	密度の異なる領域を、波が屈折しながら通過するシミュレーションです。
力学アルゴリズムによる波のシミュレーション（面上の波）
	媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。
手動で波を発生させるシミュレーション
	スライダーをマウスで動かす事により、波を発生させるシミュレーションです。
力学アルゴリズムによる波のシミュレーション（線上の波）
	媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。
二重振り子のシミュレーション
	ラグランジュ方程式を用いた、二重振り子のシミュレーションです。
ローレンツアトラクタ（GUI版）
	4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。

この階層の目次

VCSSL > 日本語トップ > コードアーカイブ > シミュレーション

[ 前へ | 目次 | 次へ ]

お知らせ

各ソフトウェアをアップデート、リニアングラフのコマンド拡張やVCSSLの英語対応など
2024/02/05 - 各ソフトの一斉アップデートの内容をお知らせします。今回は、リニアングラフのコマンド機能を大幅拡張したのがメインです。また、VCSSLのメッセージ類の英語対応も行いました。

Vnano の Ver.1.1 で実装した反復計算高速化の内側
2024/01/17 - 前回のお知らせ記事の続編です。スクリプトエンジン Vnano の Ver.1.1 において実施した高速化を、エンジン内部の実装面から掘り下げて解説します。

スクリプトエンジン Vnano の Ver.1.1 を公開、同じ内容の反復実行が劇的に高速化
2023/12/22 - スクリプトエンジン Vnano の最新版を公開しました。同じ計算式やスクリプトを反復実行する速さが、大幅に向上しました。次期版リニアングラフ3D内での活用例を交えつつ、詳しく解説します。

» お知らせ一覧

新着

[公式ガイドサンプル] 各種の図形や画像を描画する
	「VCSSL 2DCG開発ガイド」内のサンプルコードです。図形や画像などを描画します。 2022-12-16
角度の「度」とラジアンとを相互変換し、図示もするツール
	45度などの「度」の値と、ラジアンの値とを相互に変換できるツールです。対応する角度の図示もできます。 2022-11-22
[公式ガイドサンプル] 立体モデルを回転させるアニメーション
	「VCSSL 3DCG開発ガイド」内のサンプルコードです。立体モデルを、アニメーションで回転させ続けます 2022-11-18
[公式ガイドサンプル] 立体モデルを回転させる
	「VCSSL 3DCG開発ガイド」内のサンプルコードです。立体モデルを、指定した角度だけ回転させます。 2022-11-17
[公式ガイドサンプル] 四則演算の計算順序や、カッコを使った順序指定
	「VCSSLスタートアップガイド」内のサンプルコードです。四則演算が計算される順序や、カッコを使って順序を指定する方法を例示します。 2022-10-31

Twitterでも
お知らせを配信しています！

プログラミング言語VCSSLやVnano、およびこのウェブサイトは、RINEARNが開発/運営しています。

プログラミング言語VCSSL 公式サイト

このウェブサイトでは、広告配信とトラフィックデータ解析にCookieを使用しています。このサイトをご利用頂くにはプライバシーポリシーをご確認ください。
このサイトはRINEARNが管理・運営を行っています。サイト内に掲載されている文章・画像の著作権は、原則としてRINEARNおよびその代表に帰属します。

二重振り子のシミュレーション

スポンサーリンク

使用方法

ダウンロードと展開（解凍）

プログラムの起動

Windows をご使用の場合

Linux 等をご使用の場合

操作方法

題材解説

二重振り子のモデル

ニュートン方程式で素直に扱うと難しいので、解析力学で扱う

ラグランジアンとは

運動エネルギー T を求める

位置エネルギー U を求める

ラグランジアンとラグランジュ方程式

行列化

逆行列を用いて、行列式を解く

オイラー法

スポンサーリンク

コード解説

コーディング

メインループ型の処理フロー

テンプレート

ヘッダ領域の記述

グローバル変数の宣言

中身が空の関数は、目的に応じた処理を独自に実装する

init関数 ― 初期化処理の実装

paint関数 ― 描画処理の実装

update関数 ― 更新処理の実装

ライセンス

この記事中の商標などについて