» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら
ローレンツアトラクタ（ファイル出力版）
このプログラムは、ローレンツ方程式を4次ルンゲ=クッタ（RK4）法で数値的に解き、 その解曲線を描くVCSSLプログラムです。 描かれる解曲線は、「 ローレンツアトラクタ 」として非常に有名です。
なお、このプログラムには、以下のGUI版も存在します。ローレンツアトラクタについて詳しい内容については、 GUI版のほうをご参照下さい。
上のGUI版では、ユーザーの操作に従って、座標値をリアルタイムで計算しながら描画させていました。 一般にローレンツアトラクタを描いて楽しむ分には、そのほうが便利です。 それに対して今回のプログラムでは、座標値を一度ファイルに書き出してから、 グラフソフトでそのファイルをプロットさせるという、数値計算において一般的な、シンプルなものとなっています。
計算結果をファイルに書き出す事で、VCSSL付属のグラフ機能だけではなく、 一般のグラフソフトで解析する事もできますし、出力結果をスクリプトなどで加工する事も容易になり、 より汎用的になります。 反面、手軽にパラメータを変えて楽しむような場合にはやや面倒になります。 そういった場合はやはりリアルタイムで計算と可視化を行う GUI版 のほうが便利です。
スポンサーリンク
使用方法
ダウンロードと展開（解凍）
まず、PC（スマホは未対応）で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。
その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開（解凍）してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。
» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…
プログラムの起動
Windows をご使用の場合
上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください：
もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。
正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。
Linux 等をご使用の場合
ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください：
（プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行）
» javaコマンドが使用できない等のエラーが表示される場合は…
操作方法
3Dグラフの画面上で、下記のマウス操作を行えます。
- 左ドラッグ … 視点の回転
- 右ドラッグ … 視点の平行移動
- ホイールスクロール … 拡大/縮小
パラメータの変更
波に関する各種パラメータは、プログラムの先頭領域にまとめて記載されています。 それらの値を書き換えると、色々と条件を変えてシミュレーションを楽しめます。
題材解説
ローレンツアトラクタについての詳しい解説は、下記のGUI版をご参照下さい。
コード解説
コード全体
このプログラムのコード全体は、以下の通りです。
coding UTF-8;
import Math;
import tool.Graph3D;
const string FILE_NAME = "lorenz.dat3d"; // 出力ファイル名
const double S = 10.0; // ローレンツ方程式のパラメータ: σ
const double B = 4.7; // ローレンツ方程式のパラメータ: β
const double R = 100.7; // ローレンツ方程式のパラメータ: ρ
const double TIME_MAX = 8.0; // 終了時刻
const double DT = 1.0E-4; // 1サイクルの時間進行
const int GRAPH_POINTS = 2400; // グラフプロット点数
const int LOOPS = TIME_MAX / DT; // 必要ループ回数
void main(){
// 軌道の座標値
double x = 0.1;
double y = 0.1;
double z = 0.1;
// RK4法で必要な変数
double k1_x, k1_y, k1_z;
double k2_x, k2_y, k2_z;
double k3_x, k3_y, k3_z;
double k4_x, k4_y, k4_z;
// TSV書き込みモードでファイル「map.dat2d」を開く
int fileID = open(FILE_NAME, "wtsv");
println("計算を開始します...");
// RK4（4次ルンゲ=クッタ）法で座標値を計算
for( int i=0; i<LOOPS; i++ ){
k1_x = fX(x, y, z, S, B, R);
k1_y = fY(x, y, z, S, B, R);
k1_z = fZ(x, y, z, S, B, R);
k2_x = fX(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k2_y = fY(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k2_z = fZ(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k3_x = fX(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k3_y = fY(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k3_z = fZ(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k4_x = fX(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
k4_y = fY(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
k4_z = fZ(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
x += (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x)*DT/6.0;
y += (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y)*DT/6.0;
z += (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z)*DT/6.0;
// ファイルに座標値を書き出し
if( i % (LOOPS/GRAPH_POINTS) == 0 ){
writeln(fileID, x, y, z);
}
}
// ファイルを閉じる
close(fileID);
println("計算が終了しました。グラフにプロットします。");
// ファイルをグラフにプロット
int graph = newGraph3D();
setGraph3DFile(graph, FILE_NAME);
}
/**
* X方向の変化率を与える関数
*/
double fX(double x, double y, double z, double s, double b, double r){
return -s*x + s*y;
}
/**
* Y方向の変化率を与える関数
*/
double fY(double x, double y, double z, double s, double b, double r){
return -x*z + r*x - y;
}
/**
* Z方向の変化率を与える関数
*/
double fZ(double x, double y, double z, double s, double b, double r){
return x*y - b*z;
}
このプログラムのコードは上記の通り、典型的な数値計算プログラムの内容となっています。
以下では、このコードの各部を順を追って解説します。
ヘッダ部 - ライブラリのインポートなど
まずコード先頭のヘッダ部では、文字コードの明示（任意）と、必要なライブラリのインポートなどを行っています。
import Math;
import tool.Graph3D;
「 Math 」は各種数学関数を提供する標準ライブラリで、 「 tool.Graph3D 」は3次元グラフ描画機能を提供するライブラリです。
グローバル領域 - 計算のパラメータなど
続いてグローバル領域ですが、ここではローレンツ方程式を数値的に解く際のパラメータなどを定数で定義しています。
const string FILE_NAME = "lorenz.dat3d"; // 出力ファイル名
const double S = 10.0; // ローレンツ方程式のパラメータ: σ
const double B = 4.7; // ローレンツ方程式のパラメータ: β
const double R = 100.7; // ローレンツ方程式のパラメータ: ρ
const double TIME_MAX = 8.0; // 終了時刻
const double DT = 1.0E-4; // 1サイクルの時間進行
const int GRAPH_POINTS = 2400; // グラフプロット点数
const int LOOPS = TIME_MAX / DT; // 必要ループ回数
実はこれらは全てmain関数のローカルスコープ内に入れてしまっても動くので、 ある意味で無駄にグローバル領域を使っています。
しかし GUI版 や、 モジュール外から初期値を操作して何度も再計算させるような場合では、 パラメータはグローバル領域に置いておく必要があります。 今回は短い単純な計算コードで、あまり色々気をつけて考えても大した利点は無いので、 ここではグローバルに置きっぱなしにしています。
S, B, R がローレンツ方程式のパラメータです。なお、ローレンツ方程式は以下のような形をした1階の常微分方程式です。
\[ \frac{dx}{dt} = -sx + sy \] \[ \frac{dy}{dt} = -xz + rz - y \] \[ \frac{dy}{dt} = xy - bz \]
今回のプログラムでは、この微分方程式の解曲線を RK4（ 4次ルンゲ=クッタ法 ） で求めています。RK4は科学技術計算で実際に日常的に使われている高精度な手法で、 よく知られている オイラー法 などの強力版に相当するものです。
RK4は離散化幅 DT （ DTやdtと書くのは時間微分の方程式の場合で、一般の場合には h や Δ と書かれる事も多い ） に対して4次の精度を誇るアルゴリズムで、つまりは DT を1桁細かくすると、結果の精度は 1万倍向上します。 従って、絶対的な精度が要求される場合に威力を発揮しますが、 それでありながら1ループの計算コストもそこそこに納まるため、便利で扱いやすいアルゴリズムです。
パラメータ定数の TIME_MAX は終了時刻、DT は計算の際における時間の離散化幅です。 GRAPH_POINTS は3Dグラフに何点プロットする点数で、つまりは座標値ファイルに書き出す点数です。 あまり細かすぎても重い上に潰れて意味が無いため、ある程度省いて書き出すようにしています。
計算部
それでは中心である計算部を見ていきましょう。今回は単純なので、main関数の中にそのまま書いています。
// 軌道の座標値
double x = 0.1;
double y = 0.1;
double z = 0.1;
// RK4法で必要な変数
double k1_x, k1_y, k1_z;
double k2_x, k2_y, k2_z;
double k3_x, k3_y, k3_z;
double k4_x, k4_y, k4_z;
// TSV書き込みモードでファイル「map.dat2d」を開く
int fileID = open(FILE_NAME, "wtsv");
println("計算を開始します...");
// RK4（4次ルンゲ=クッタ）法で座標値を計算
for( int i=0; i<LOOPS; i++ ){
k1_x = fX(x, y, z, S, B, R);
k1_y = fY(x, y, z, S, B, R);
k1_z = fZ(x, y, z, S, B, R);
k2_x = fX(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k2_y = fY(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k2_z = fZ(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
k3_x = fX(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k3_y = fY(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k3_z = fZ(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
k4_x = fX(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
k4_y = fY(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
k4_z = fZ(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
x += (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x)*DT/6.0;
y += (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y)*DT/6.0;
z += (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z)*DT/6.0;
// ファイルに座標値を書き出し
if( i % (LOOPS/GRAPH_POINTS) == 0 ){
writeln(fileID, x, y, z);
}
}
// ファイルを閉じる
close(fileID);
前半部では、まず座標変数(x,y,z)とRK4に必要な変数を用意し、open 関数でファイルを開いています。 その後は、forループの中でRK4のアルゴリズムを回して、一歩ずつ端から解曲線の座標を求め、 writeln 関数でファイルに書き出しています。
forループ内の処理は典型的なRK4のコードで、何も特殊な事は行っていません。 fX, fY, fZはそれぞれ \(dx/dt\), \(dy/dt\), \(dz/dt\) の値、つまりはローレンツ方程式の右辺を返す関数で、 その実装は以下の通りです：
/**
* X方向の変化率を与える関数
*/
double fX(double x, double y, double z, double s, double b, double r){
return -s*x + s*y;
}
/**
* Y方向の変化率を与える関数
*/
double fY(double x, double y, double z, double s, double b, double r){
return -x*z + r*x - y;
}
/**
* Z方向の変化率を与える関数
*/
double fZ(double x, double y, double z, double s, double b, double r){
return x*y - b*z;
}
この変化率関数を用いて、現在の座標から次の時刻（DT後）の座標を求め、 それがforループの1サイクルになっています。次々とDTずつ次の時刻の座標が求まり、 最終的に解曲線が求まるわけです。 この流れは基本的にオイラー法などと変わりません。 異なるのは、1ループの中で少しずつ半端な場所を取りながら何度も導関数 fX, fY, fZ を“見て”いるのと、最終的にそれらに重みをつけて平均している点です。
グラフプロット部
最後に、書き出したファイルをグラフにプロットする部分です。
// ファイルをグラフにプロット
int graph = newGraph3D();
setGraph3DFile(graph, FILE_NAME);}
newGraph3D 関数は、3Dグラフソフトを起動し、 そのID（識別番号）を整数で返します。setGraph3DFile 関数は、 そうして起動したグラフソフトに、座標値ファイルを開いてプロットさせる関数です。
GUI版 では、setGraph3DData 関数を使って、 ファイルでは無く、全座標値を格納する float 配列をグラフにプロットさせていました。 速度的にはそちらのほうが遥かに速いのですが、今回のようにファイルを経由する事で、 他のグラフソフトや解析ツールを併用する事などができるようになります。 また、今回は一瞬で計算が終わりますが、長い時間がかかるような場合にも、 毎回計算していては大変なので、座標値をファイルに書き出して使う事になります。 そういった場合は、今回のように一旦ファイルに書き出して、setGraph3DFile 関数でプロットさせるのが便利です。
ライセンス
この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 氏の米国およびその他の国における商標または登録商標です。
- その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。
|
Vnano版 | ローレンツ方程式を数値的に解くプログラム |
|
 |
ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。 |
|
波の干渉（面上の円形波）のアニメーション表示 |
|
 |
面上の円形波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
|
円形波のアニメーション表示 |
|
 |
振幅・波長・周期をスライダ―で操作しながら、円形波のグラフをアニメーションで見られるプログラムです。 |
|
波の干渉（線上の正弦波）のアニメーション表示 |
|
 |
線上（1次元の）の正弦波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
|
正弦波のアニメーション表示 |
|
 |
振幅・波長・周期をスライダ―で操作しながら、正弦波のグラフをアニメーションで見られるプログラムです。 |
|
凹レンズを通過する波のシミュレーション |
|
 |
凹レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
|
凸レンズを通過する波のシミュレーション |
|
 |
凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
|
乱雑な密度分布における波のシミュレーション |
|
 |
密度分布が乱雑な媒質中における、波の伝播のシミュレーションです。 |
|
ローレンツアトラクタ（ファイル出力版） |
|
 |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |
|
波の屈折のシミュレーション |
|
 |
密度の異なる領域を、波が屈折しながら通過するシミュレーションです。 |
|
力学アルゴリズムによる波のシミュレーション（面上の波） |
|
 |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
|
手動で波を発生させるシミュレーション |
|
 |
スライダーをマウスで動かす事により、波を発生させるシミュレーションです。 |
|
力学アルゴリズムによる波のシミュレーション（線上の波） |
|
 |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
|
二重振り子のシミュレーション |
|
 |
ラグランジュ方程式を用いた、二重振り子のシミュレーションです。 |
|
ローレンツアトラクタ（GUI版） |
|
 |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |
