ローレンツアトラクタ（ファイル出力版）

このプログラムは、ローレンツ方程式を4次ルンゲ=クッタ（RK4）法で数値的に解き、 その解曲線を描くVCSSLプログラムです。 描かれる解曲線は、「 ローレンツアトラクタ 」として非常に有名です。

なお、このプログラムには、以下のGUI版も存在します。ローレンツアトラクタについて詳しい内容については、 GUI版のほうをご参照下さい。

• ローレンツアトラクタ（GUI版）

上のGUI版では、ユーザーの操作に従って、座標値をリアルタイムで計算しながら描画させていました。 一般にローレンツアトラクタを描いて楽しむ分には、そのほうが便利です。 それに対して今回のプログラムでは、座標値を一度ファイルに書き出してから、 グラフソフトでそのファイルをプロットさせるという、数値計算において一般的な、シンプルなものとなっています。

計算結果をファイルに書き出す事で、VCSSL付属のグラフ機能だけではなく、 一般のグラフソフトで解析する事もできますし、出力結果をスクリプトなどで加工する事も容易になり、 より汎用的になります。 反面、手軽にパラメータを変えて楽しむような場合にはやや面倒になります。 そういった場合はやはりリアルタイムで計算と可視化を行う GUI版 のほうが便利です。

スポンサーリンク

---

使用方法

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

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

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

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

プログラムの起動

Windows をご使用の場合

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

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

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

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

Linux 等をご使用の場合

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

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

操作方法

3Dグラフの画面上で、下記のマウス操作を行えます。

• 左ドラッグ … 視点の回転
• 右ドラッグ … 視点の平行移動
• ホイールスクロール … 拡大/縮小

パラメータの変更

波に関する各種パラメータは、プログラムの先頭領域にまとめて記載されています。 それらの値を書き換えると、色々と条件を変えてシミュレーションを楽しめます。

題材解説

ローレンツアトラクタについての詳しい解説は、下記のGUI版をご参照下さい。

• ローレンツアトラクタ（GUI版）

コード解説

コード全体

このプログラムのコード全体は、以下の通りです。

このプログラムのコードは上記の通り、典型的な数値計算プログラムの内容となっています。

以下では、このコードの各部を順を追って解説します。

ヘッダ部 - ライブラリのインポートなど

まずコード先頭のヘッダ部では、文字コードの明示（任意）と、必要なライブラリのインポートなどを行っています。

「 Math 」は各種数学関数を提供する標準ライブラリで、 「 tool.Graph3D 」は3次元グラフ描画機能を提供するライブラリです。

グローバル領域 - 計算のパラメータなど

続いてグローバル領域ですが、ここではローレンツ方程式を数値的に解く際のパラメータなどを定数で定義しています。

実はこれらは全て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関数の中にそのまま書いています。

前半部では、まず座標変数(x,y,z)とRK4に必要な変数を用意し、open 関数でファイルを開いています。 その後は、forループの中でRK4のアルゴリズムを回して、一歩ずつ端から解曲線の座標を求め、 writeln 関数でファイルに書き出しています。

forループ内の処理は典型的なRK4のコードで、何も特殊な事は行っていません。 fX, fY, fZはそれぞれ \(dx/dt\), \(dy/dt\), \(dz/dt\) の値、つまりはローレンツ方程式の右辺を返す関数で、 その実装は以下の通りです：

この変化率関数を用いて、現在の座標から次の時刻（DT後）の座標を求め、 それがforループの1サイクルになっています。次々とDTずつ次の時刻の座標が求まり、 最終的に解曲線が求まるわけです。 この流れは基本的にオイラー法などと変わりません。 異なるのは、1ループの中で少しずつ半端な場所を取りながら何度も導関数 fX, fY, fZ を“見て”いるのと、最終的にそれらに重みをつけて平均している点です。

グラフプロット部

最後に、書き出したファイルをグラフにプロットする部分です。

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 氏の米国およびその他の国における商標または登録商標です。
• その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。