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

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

使用方法

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

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

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

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

プログラムの起動

Windows をご使用の場合

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

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

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

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

Linux 等をご使用の場合

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

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

操作方法

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

題材解説

二重振り子のモデル

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

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

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

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

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

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

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

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

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

ラグランジアンとは

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

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

\[ \mathcal{L} = T - U \]

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

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

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

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

\[ V_{1}^2 = V_{1,x}^2 + V_{1,y}^2 = (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}^2 = V_{2,x}^2 + V_{2,y}^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つめの振り子の支点の高さにとりましょう。 そして、支点から鉛直上向きにy軸をとります。 すると、ボールが支点よりも下にある時に、位置エネルギーがマイナスになりますが、特に問題はありません。

この場合、このモデルの位置エネルギー U は：

\[ 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\) は重力加速度です。なお、上の第1項が 1 つめのおもりの位置エネルギー、第2項が 2 つめのおもりの位置エネルギーです。

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

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

\[ \mathcal{L} = T - U = \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 \mathcal{L}}{\partial\dot\theta_{1}}-\frac{\partial \mathcal{L}}{\partial\theta_{1}}=0 \] \[ \frac{d}{dt}\frac{\partial \mathcal{L}}{\partial\dot\theta_{2}}-\frac{\partial L}{\partial\theta_{2}}=0 \]

分野によっては、\( \partial \mathcal{L} / \partial\theta \) といった表記は見慣れないものかもしれません。 これは偏微分と呼ばれるもので、\( \theta \) 以外の変数をすべて定数と見なす、やや強引な微分を意味します。 \( \partial \mathcal{L} / \partial\theta \) の計算では、\( \dot\theta \) や \( \ddot\theta \) も定数と見なします。 つまり、特に何も考えず、\( \theta \) だけに着目して微分すれば良いのです（普通の微分では、もっと変数間の依存関係を考えないといけません）。

同様に \( \partial \mathcal{L} / \partial\dot\theta \) では、今度は \( \dot\theta \) のみを変数と見なし、\( \theta \) や \( \ddot\theta \) などはすべて定数と見なして微分します。

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

\[ (M_{1} + M_{2})L_{1}^{2}\ddot\theta_{1} + M_{2}L_{1}L_{2} \ddot\theta_{2} \cos(\theta_{1}-\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} \ddot\theta_{1} \cos(\theta_{1}-\theta_{2}) + 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) を導けました。

「導けました」って、これで一体どうシミュレーションするの…？ というのはこの後で説明していきます。

行列化

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

\[ \ddot\theta_{1}=... 　　(3) \] \[ \ddot\theta_{2}=... 　　(4) \]

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

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

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

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

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

\( \bigg( \) \( A_{1} \) \( B \) \( \bigg )\) \( B \) \( A_{2} \)

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

\( = \)

\( \bigg( \) \( D_1 \) \( \bigg) \) \( D_2 \)

(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_{2}=M_{2}L_{1}L_{2}\dot\theta_{1}\sin(\theta_{1}-\theta_{2})-M_{2}gL_{2}\sin\theta_{2} 　　(10) \]

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

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

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

\( \bigg( \) \( a \) \( b \) \( \bigg)^{-1} \) \( c \) \( d \)

\( = \)

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

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

\( \bigg( \) \( A_{1} \) \( B \) \( \bigg)^{-1} \) \( B \) \( A_{2} \)

\( = \)

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

と求まります。

ここで扱っている 2×2 の場合だけでなく、3×3 の行列についても逆行列を求める公式が存在します。4×4 でも、だいぶ複雑になりますが存在はします。 さらに、コンピューターで数値的に逆行列を求めるアルゴリズムも完成されており、行列演算ライブラリに当然のように入っています。 次元数が多い場合は、そうやって数値的に求めるのも現実的な手でしょう。 ただしその場合、上の A, B, D の値は \( \theta_i \) や \( \dot\theta_i \)に依存しているため、行列は瞬間瞬間で変化する（逐次的に更新が必要になる）事を忘れないよう注意です。

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

\( \bigg( \) \( A_{1} \) \( B \) \( \bigg)^{-1} \) \( B \) \( A_{2} \)

\( \bigg( \) \( A_{1} \) \( B \) \( \bigg) \) \( B \) \( A_{2} \)

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

\( = \)

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

\( \bigg( \) \( D_{1} \) \( \bigg) \) \( D_{2} \)

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

\( \bigg( \) \( 1 \) \( 0 \) \( \bigg) \) \( 0 \) \( 1 \)

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

\( = \)

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

よって解は：

\[ \ddot\theta_{1}=\frac{A_{2}D_{1}-BD_{2}}{A_{1}A_{2}-B^{2}} \] \[ \ddot\theta_{2}=\frac{A_{1}D_{2}-BD_{1}}{A_{1}A_{2}-B^{2}} \]

と、このように求まりました。なお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) \]

ここで「\(\approx\)」は近似的に近い関係（厳密なイコールではない）を表します。

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

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

\[ \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の値も \(\theta\) や \( \dot\theta \) に依存しているため、これらの値も毎時刻、(6)〜(10)式を用いて再計算してやる必要があります。

コード解説

コーディング

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

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

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

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

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

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

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

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

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

テンプレート

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

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

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

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

ヘッダ領域の記述

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

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

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

グローバル変数の宣言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

paint関数 ― 描画処理の実装

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

ここで drawPoint 関数などに渡すのは画面座標系の値であり、これは画面の上端を原点として、右向きにX軸、下向きにY軸をとったものです。単位はピクセルです。

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

update関数 ― 更新処理の実装

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

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

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

\[ \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) \]

ここで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_{2}=M_{2}L_{1}L_{2}\dot\theta_{1}\sin(\theta_{1}-\theta_{2})-M_{2}gL_{2}\sin\theta_{2} 　　(10) \]

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

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

ライセンス

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

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

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

参考文献

[1]: ランダウ＝リフシッツ理論物理学教程 力学（増訂第3版）, 東京図書株式会社, 1974年

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

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