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

Vnano版 | ローレンツ方程式を数値的に解くプログラム

ローレンツ方程式( » 詳しくは 題材解説 参照 )を4次のルンゲ=クッタ法によって数値的に解き、求まった解曲線を、グラフ描画用のデータとして出力するプログラムです。 グラフの描画方法についても解説します。 描かれる解曲線は、いわゆるローレンツアトラクタとして有名です。

コードは Vnano (下記参照) で記述していますが、コード解説 の項目において、C言語/C++/Java言語で書いたコードも配布しています。 また、Vnano は VCSSL のサブセットであるため、拡張子を変えればそのままVCSSLのコードとしても動きます。

Vnano とは?

Vnano (VCSSL nano) は、Java®言語製のソフトウェア内に組み込んでスクリプト処理機能として使える、小型スクリプトエンジン&言語です。 C言語系のシンプルな文法を採用し、計算/解析ソフト等での使用も想定して数値演算速度を重視しています。 例えば関数電卓の「 RINPn 」では、Vnano の計算コードを実行したり、関数定義に使う事ができます。

スポンサーリンク


使用方法

※ このコードの記述言語である Vnano は、現在はまだ開発後半のβ版の段階です。あらかじめご了承ください。 このコードの実行における動作検証は済んでいますが、コード内容を大々的に書き換える場合などには、上記の点に留意が必要です。

ダウンロードと展開(解凍)

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

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

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

※ Linux 等でも右クリックメニューから展開できますが、ファイル名が文字化けする場合、コマンド端末でZIPファイルのある場所まで cd した上で「 unzip -O cp932 ZIPファイル名 」で展開してみてください。

プログラムの起動

Microsoft® Windows® をご使用の場合

上記の通りにZIPファイルを展開したフォルダ内にある、 「 VCSSL.bat(種類はバッチファイル) 」をダブルクリック実行してください。 もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor.bat 」を実行してください。

※ もし、ここでセキュリティ警告メッセージが出て起動できない場合は、 ダウンロードと展開 の項目に書かれている手順通りに、ZIPファイルの展開をやり直してみてください。

正常に起動できると、最初にメモリー使用量や、(必要な場合のみ)Java®実行環境を自動で入手するか 等を尋ねられるので、適時答えると、プログラムが起動します。2回目以降はすぐに起動します。

※ ここで入手したJava®実行環境は、ZIPファイルを展開した中の「 jre 」フォルダ内にダウンロードされ、このプログラムの実行のみに使用されます。PC全体に影響する形でインストールされる事はありません。

Linux 等やその他のOSをご使用の場合

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

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

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

起動後

グラフ化用データの表示、またはファイルへの保存

起動すると、まず黒いウィンドウが表示されて計算処理が始まり、 その画面上にグラフ用データが書き出されていきます:

データが表示される画面の様子

データ内には、グラフに描く座標値が、一行につき一点ずつ「 x    y    z 」とタブ空白区切りで記載されています。

計算が完了すると、黒い画面へのデータ書き出しも止まります(標準では一瞬で終わります)。 ここで 画面上部のメニューバーから「 File 」>「 Save File ... 」を選択すると、 表示されたデータをファイルに保存できます。

表示/保存したデータからグラフを描く方法については、グラフの描画方法 の項目をご参照ください。

グラフの描画方法

VCSSLランタイムでの場合(今の手順での場合はこちら)

このページからダウンロードしたパッケージでは、実行用にVCSSLランタイム(※コード同梱配布用のコンパクト版)というソフトが使われます。 VCSSLランタイムでは、黒い画面に表示されたデータをグラフ化する機能が標準搭載されているため、何も用意せにグラフを描画できます。

具体的には、計算完了後、黒い画面上部のメニューバーから「 File 」>「 Plot Graph3D ... 」を選択すると…

グラフをプロットするためのメニュー項目の様子

すると、3Dグラフの画面が表示されます。標準では、グラフは点のみで描画されますが、 グラフ画面上部のメニューバーから「 オプション 」>「 線プロット 」をONにすると、 以下のように点が線でつながれたグラフが得られます:

表示されるグラフ

※ X/Y/Z軸範囲や目盛り数値の書式設定などは、グラフ画面上部のメニューバーから行えます。 グラフを資料で使う場合などには、そういった見栄えの調整を行ってください。 標準では上のように、単純にX/Y/Z軸範囲はデータの最大最小値そのものに設定されます。 これはデータのありのままを把握するのには好都合なので、自動で見栄え用の調整等はされないようになっています。

このグラフが、プログラム内で指定されたパラメータ・初期値の条件下で、4次のルンゲ=クッタ法によって数値的に解いて求めた、ローレンツ方程式の解曲線をプロットしたものとなっています。

ここでグラフ画面上部のメニューバーから、「ツール」 > 「アニメーション」を選択し、「START」ボタンを押すと、 解曲線を始点(初期値)から終点に向かって描いていくようなアニメーションを行えます。

実際のアニメーション例はこちら(約6.5MB)です。

関数電卓 RINPn を使っている場合

関数電卓ソフトの RINPn では、Vnanoのプログラムは実行できますが、グラフ表示機能は現時点では標準搭載されていません。 そこで、別途グラフソフトを用意し、組み合わせて使用します。 今回のプログラムが出力するデータは、結構スタンダードな形式なので、gnuplot などの色々なグラフソフトで描画できます

ここでは例として、VCSSL や Vnano、RINPn と同じ開発元(つまりうちです)のグラフソフトである「 リニアングラフ 」を使ってみましょう。 2D版と3D版がありますが、今回は3D版を使用します:

リニアングラフ3D   公式サイト(日本語ページ)
https://www.rinearn.com/ja-jp/graph3d/

リニアングラフを起動すると、先に扱ったVCSSLランタイムの場合と全く同じグラフ画面が立ち上がります。 それもそのはず、VCSSLランタイムに標準搭載されていたグラフ機能の中身もリニアングラフだからです。

グラフの描画方法は簡単です。 まずは普通に RINPn で今回のプログラムを実行しましょう。 プログラムのファイル「 LorenzAttractor.vnano 」を RINPn のフォルダ内に置いて、 RINPn の電卓画面のINPUT項目に、このファイル名をそのまま入力します。 すると黒い画面が表示され、そこにグラフ用のデータが出力されます:

RINPnでデータが出力される様子
※ 画面の様子はOSや環境によって多少異なります。 数値の桁数がところどころ不揃いなのは、プログラム側で丸めを行っていないためで、VCSSLランタイムの項目での同様の説明をご参照ください。

この黒い画面のデータを全行選択(マウスでドラッグしてもいいですが、Ctrl+Aキーの同時押しもで可能です)した状態で、右クリックメニューからコピーし、 そしてリニアングラフ3Dの画面上で右クリックして「 データの貼り付け 」を行うと、グラフが表示されます:

表示されるグラフ

VCSSLランタイムの場合と同じグラフですね。 どちらも内部で同じスクリプトエンジン(Vnanoエンジン)を使って計算処理を行っていて、 それを実質同じグラフソフトで描画しているため、どちらを使っても完全に同一のグラフが得られます。

RINPn を使ってコマンドラインで実行する場合

RINPn は、「 bin 」フォルダのパスをOSの環境変数 Path/PATH に登録すると「 rinpn 」コマンドが使用できます。 その場合はリニアングラフの bin フォルダのパスも一緒に登録しておくと、こちらも「 ring3d 」コマンド(2Dの場合は ring2d )が使えるようになり、 両者を組み合わせるのに便利です。実際、今の場合の例では、コマンドラインで:

cd "LorenzAttractor.vnanoがある場所"
rinpn LorenzAttractor.vnano > data.txt
ring3d data.txt

とすると、計算結果がファイル data.txt に保存された上で、その内容がグラフとして表示されます。 表示されるグラフの内容は、上で掲載したものと全く同じです。

RINPn をコマンドラインで使いつつ、描画に gnuplot を使いたい場合

余談ですが、グラフ描画に gnuplot を使いたいという場合は:

cd "LorenzAttractor.vnanoがある場所"
rinpn LorenzAttractor.vnano > data.txt
gnuplot
(ここでgnuplotが起動します)
set xlabel "X"
set ylabel "Y"
set zlabel "Z"
splot "data.txt" with lp palette pt 7 ps 0.5
(使い終わったら)
quit

などとします。「splot 〜」の行がメインのプロットコマンドです。「 with lp 」が線とマーカー(点)でのプロット、「 palette 」がグラデーション彩色する指定、「 pt 7 」がマーカーを丸い点にする指定(環境やバージョンによっては別の数字かもしれません)、「 ps 0.5 」がそのサイズ指定です。 表示されるグラフは:

表示されるグラフ

この通りです。細かい見栄えは違いますが、グラフの形としては先ほどのリニアングラフ3Dの場合と同じですね。 このように、今回のコードで出力されるような書式のデータ(各行に1座標点ずつ x   y   z と並べたデータ)は、色々なグラフソフトでプロットできます。

なお、gnuplot はグラフソフトの大御所という感じのソフトで、計算関連ではこのソフトに関わらずに過ごし続けるのがむしろ難しいようなソフトです。 印刷時に綺麗な形式で画像を出力できるので、そういう図を作る場面などでは特に重宝します。 慣れないとコマンド指定が難しく感じるかもしれませんが、コマンド列をファイルに書いたものを load コマンドで読み込めるので、作業開始時に一度書けば使い回せます。 各種OS用のものがあるので、計算やグラフ描画が面白いと感じた方は、PCに入れておいて損はないと思います。

スポンサーリンク


題材解説

ローレンツ方程式

今回のプログラムは、「 ローレンツ方程式 」を数値的に解いて、その解の曲線をグラフにプロットするためのものです。 そこで、プログラムのコード内容の解説に入る前に、今回のテーマであるローレンツ方程式について触れておきましょう。

ローレンツ方程式は、アメリカのMITの気象学者であった エドワード・ローレンツ 氏が1963年(提出は1962年)に出した この論文(アメリカ気象学会のサイト内で一般公開されています) に初登場する、以下の形の微分方程式です(係数の文字は r→\(\rho\)、b→\(\beta\) に変えています):

\[ \frac{dx}{dt} = - \sigma x + \sigma y \] \[ \frac{dy}{dt} = - x z + \rho x + y \] \[ \frac{dz}{dt} = x y - \beta z \]

一応は Wikipedia にもページがありますが、 日本語版のページ はやや内容が少なく、 英語版のページ の方が充実しています。方程式の背景にある歴史やモデル等についても掘り下げて知りたいという方は、後者のページや元論文を併せて参照してみてください。

解曲線

さて、上のローレンツ方程式は、\(x, y, z\) 座標の時間 \(t\) に関する微分の形(= 時間経過に対する瞬間的な変化)を与えている方程式です。 なので、要は3次元の点の\(x, y, z\) 座標が、ある時刻 \(t\) の瞬間に、それぞれどう変化するかという、「動き」のルールを与えているとイメージする事ができます。

つまり最初の時刻 \(t = 0\) において、どこか適当な位置(\(x, y, z\) の初期値)に点を置くと、 その点は時間経過に伴って、上の方程式に従って動いていきます。 結果、その点の「動き」は上の方程式を満たしている = 方程式の解になっており、描かれる軌道を解曲線と呼びます。

実際にローレンツ方程式のパラメータ \( \sigma, \rho, \beta \) をほどほどにいい感じの値(プログラム内で書き換えられますので、ぜひ色々変えて遊んでみてください)に設定して、 ほどほどに適当な初期値 \(x,y,z\) からスタートした際の解曲線が、先ほど見たグラフの通りになります:

表示されるグラフ

上図の通り、パラメータと初期値によっては、ローレンツ方程式の解曲線には独特の「 ∞ 」の字のようなパターンが描かれます。 この独特な構造は、アトラクタという種類の構造の一つで、一般にローレンツ・アトラクタと呼ばれています。

さらに面白いのが、点の動きです。 以下のアニメーションの通り、点はその「∞」のパターンの片方の輪っかを回っているかと思えば、 突然もう片方の輪っかに移ったり、また突然戻ってきたり…と、気まぐれに行き来します。

アニメーションの様子

この「 2つの輪っかの間の気まぐれな行き来 」は、いつまでも続くのですが、それがこの解曲線の面白い特徴の一つです。 だって、「 動きのルール 」を司っているローレンツ方程式そのものは極めてシンプルな形で、乱数なども一切入っていない(=決定論的な)わけです。 なのに、その振る舞いは、まるで内部で乱数でも振っているかのように、気まぐれに見えるわけです。面白くないですか?

初期値鋭敏性

ところで、いま、ある場所にいる点に注目していて、「 その点がしばらく後にどのあたりにいるか 」を予測したいとしましょう。 細かい精度ではなく、雑におおざっぱに、その点が次にどちらの輪っかを回って、その次にどちらの輪っかを回って、その次に… が予測できればいいとします。 そんなもの、上で行った調子でコンピュータに計算させれば、どこまでも簡単に予想できるように思えませんか?

でも、実はそう簡単な事ではありません。というのも、上の解曲線には、最初の点の位置がほんのちょっとでもずれると、その差が成長して、 やがて 2 つの輪っかの間を行き来するパターンもどんどんずれていく…という現象が見られるからです。

実際にその様子をアニメーションにしたのが以下です:

アニメーションの様子

上の動画は、初期状態で近い(しかし少しだけ離れた)位置に置いた 2 つの点の動きを、それぞれ赤色と緑色の点として、重ねてアニメーションさせているものです。

※ 今回のプログラムの出力結果を一旦テキストファイルに保存し、 リニアングラフ3Dのメニューバーから「 ファイル 」> 「 ファイルを開く 」でそのファイルを開くと、保存したデータをグラフに描画できます。 その際、複数のファイルを開くと、複数のグラフを重ねて描画できるため、パラメータを変えて計算・保存した結果を比較する事ができます。 その状態で、アニメーションメニューを使用すると、上のようなアニメーションを再生できます。

上のアニメーションの通り、 最初は非常に近くで一緒に動いていた 2 つの点が、やがて離れ、最後はそれぞれ異なるパターンで輪っかの間を行き来する ようになってしまっていますね。

このように、初期状態のわずかな差が成長して、最終的に大きな差になってしまうのは、 ローレンツ方程式の下での点の動きが、パラメータや初期値を適切に選んでやると、初期値鋭敏性という性質を示すためです。

力学系が初期値鋭敏性を示す状況(※同じ系でも安定な領域と不安定な領域が混在していたりもします)では、 このような近傍の点との間の差が、時間経過に伴って指数関数的に成長します。 指数関数というと爆発的なペースで大きくなっていきますから、 2 つの点をどれだけ近くに置いたとしても、完全に同じ点でない限り、時間の問題でやがて離れてきてしまうわけです。 そしてある時から、どちらの輪っかを回るのかすらも別れてしまい、あとはもうそれぞれバラバラに気まぐれな動きで行ったり来たりするわけです。

予測の難しさ

そして、今ここで私たちが行っているような、コンピューターによる一般的な数値計算では、 無限の桁数は扱えないため、数値をある程度の桁数(ここでは大体16桁程度が目安)で切って扱っています。 そうすると、何かしらの計算を一手行うごとに、わずかな誤差は常に発生してしまいます。 それが初期値鋭敏性によって(指数関数的な勢いで!)成長して、やがて大きな違いになってしまう… というわけで、 そういう性質を持つような系の振る舞いを、ある程度以上長い時間にわたって予測するのは難しいだろう、という事がわかります (一応、ローレンツ方程式に限って言えば、色々と努力と進展はあるようですが)。

また、初期値鋭敏性は、ローレンツ方程式で表される系だけでなく、他の力学系においても普通によくある性質です。 従って、現実の何かの振る舞いを、シミュレーションで予測しようと思った際に、 もしその系が初期値鋭敏性を持っていた場合、初期状態の設定において現実との間にわずかでもズレがあると、 長期予測はほぼ不可能になります。そして、現実の何かの状態を、完全な精度で測るのは、それこそ不可能です。 つまるところ、現実的には未来を予測する事がほぼ不可能な系というものが、この世には普通にありふれているという事になります。

このような事をもう少し深く知りたい方は、 カオス理論バタフライ効果 について掘り下げて調べてみてください。

なお、ローレンツ氏は冒頭でも述べた通り気象学者で、 元論文内でのローレンツ方程式は、大気の循環/対流を扱うためのモデル化の際に登場します。 ローレンツ方程式の解そのものが物理的な3次元の流れを表しているわけではなく、それぞれ対応する物理量があるのですが、 いずれにしても、そのモデルがここまで単純な系にもかかわらず、上で見てきたような性質を持つ事がはっきりと示された事により、 結果的に長期的な気象予測(天気予報)の本質的な困難さを知る/説明する上での、重要な知見がもたらされました。

※ 筆者はこのあたりの歴史についてあまり詳しくないのですが、Wikipedia の バタフライ効果 のページをある程度信じるとすると、まずローレンツ氏はもう少し複雑なシミュレーションにおいてこの種の現象に出くわして着目し、 その後モデルを簡素化していって、ここまでシンプルな形でも同種の性質を示し得る、という所に着地して論文にまとめたようです。 なお、類する議論自体はローレンツ氏が始まりというわけではなく、例えば3体問題や可積分系/非可積分系に関する議論などはもっと以前から行われています。

現代の気象予測技術は、単純にローレンツ方程式に直接ひもづけられる話でもないでしょうが、 それでも、このローレンツ氏の研究が、エポックメイキングな大きな一歩であったと受け止められている事は、いまも変わりません(逆に当初はあまり注目されなかったようですが…)。 そのような歴史的経緯もあってか、現代でも長期的な気象予測の難しさや、カオスによる予測困難さについて一般向けに説明される際に、 定番のようにローレンツ氏のストーリーやローレンツ・アトラクタ、およびバタフライ効果(ローレンツ氏の例え話が発端)などが取り上げて語られます。

コード解説

ここからは、プログラムのコード内容を見ていきましょう。

コード全体

このプログラムのコードはVnanoで記述されています。 VnanoはC言語系のシンプルな文法を持っているので、C系の言語に触れられた事のある方なら、 コメントを参考にしながらコード内容を比較的簡単に追う or 改造する事ができると思います。

Vnano とは?

Vnano (VCSSL nano) は、Java®言語製のソフトウェア内に組み込んでスクリプト処理機能として使える、小型スクリプトエンジン&言語です。 C言語系のシンプルな文法を採用し、計算/解析ソフト等での使用も想定して数値演算速度を重視しています。 例えば関数電卓の「 RINPn 」では、Vnano の計算コードを実行したり、関数定義に使う事ができます。

一応は、C/C++/Java言語版のコードも併せて用意していますので、必要な方は下記のリンクを右クリックしてダウンロードしてください。 解説はVnanoのコードで行いますが、他の各言語への読み替えは簡単だと思います。

まずは、コード全体を見てみましょう:

| C言語版はこちら(SJIS/UTF-8) | C++版はこちら(SJIS/UTF-8) | Java言語版はこちら(SJIS/UTF-8) |
| C言語版はこちら(SJIS/UTF-8) | C++版はこちら(SJIS/UTF-8) | Java言語版はこちら(SJIS/UTF-8) |

大体70行ちょっとの、比較的短いコードですね。以下では、各部を掘り下げてみていきます。

コード先頭

まずは先頭です:

ここでは、いわゆる文字化けによる問題を防ぐため、使っている文字コードを宣言しています。

ローレンツ方程式や、それを数値的に解く際のパラメータの定義

続いて、ローレンツ方程式に登場するパラメータ \(\sigma, \rho, \beta\) や、 それを数値的に解く(解曲線を求める)際に使う値を定義しています:

DT 」は、積分や微分方程式を数値的に解く際のいわゆる 「 刻み幅 」で、いまの場合は、後述の計算ループの一周ごとに、時刻 \(t\) がどれだけ進むかを意味します。 ここでは 1.0E-5 つまり \(10^{-5}\) に設定しています。 基本的には細かい方が精度が期待できますが、極端に細かすぎてもだめで、解く対象ごとにほどほどの領域があります。 このあたりのイメージは、以下の数値積分の解説記事が参考になるかもしれません:

上の記事は、単に cos x を積分するだけという、今回よりずっと単純な内容ですが、それだけに全ての基本になっています。 上の記事で、刻み数 N を増やして Δx を細かくする事が、ここでの DT の値を細かくする事に(イメージ的に)対応します。

続く「 TIME_MAX 」は、 ローレンツ方程式を解いていく処理を、どの時点(時刻 \(t\))でやめるかを指定します。 題材解説の項目で述べた通り、このプログラムで求めるのは、 「 ある地点(初期値)に置いた点が、ローレンツ方程式に従って動いていく軌道(= ローレンツ方程式の解曲線) 」 です。 ローレンツ方程式そのものには「 終わりの時刻 」などは特に規定されていないため、どこかで処理を打ち切らないと、 プログラムはどこまでも点の動きを追いかけ続けて、求まる解曲線もどこまでも長く伸びていきます。 なので、見たい性質がほどほどにわかるくらいの時刻で処理を打ち切ります。

なお、今回のコードでは、点の動きを追いかける計算のループ回数 LOOPS を、単純に「 (int)(TIME_MAX / DT) 」として求めています。 この際の浮動小数点数の微小な演算誤差や、データ出力処理を簡素化(後述)している都合で、 グラフの最後の点が、必ずしも終了時刻の点と一致する事は保証されません。 TIME_MAX の値は、「だいたいそれくらいで計算を打ち切る」といった目安程度に考えてください。 もっと厳密に終了制御を行いたい場合は、自由にコードを改造して実装してみてください(意外と複雑になります)。

グラフ用のデータを出力するためのパラメータの定義

その後には、グラフに描画するためのデータを出力する際に、点の密度などを調整するためのパラメータが定義されています:

今回のプログラムでは、計算のループを1周するごとに時刻をDTずつ進めて、 いわば尺取虫(しゃくとり虫)のように点の動きを追いかけていきます。

その際、ループを1周するごとに点の座標値を出力してしまうと、DTを細かく刻んだ際や、TIME_MAX(終了時刻)を長くした際に、グラフの点数が多くなりすぎます。 今回は、解曲線の形や動きが分かれば十分なので、必要以上に点数が多くなると、描画の重さやデータサイズの点で不便です。

従って上では、グラフに大体何点くらい描かれるようにしたいかを 「 OUTPUT_POINTS 」として定義し、 そこから逆算して、計算の何ループごとに値を出力するかを OUTPUT_INTERVAL として求めています。

ただ、上の計算では、LOOPS が OUTPUT_POINTS で割り切れない場合には誤差が生じますし、 そもそも先述した LOOPS の計算自体も、ある程度の誤差を見込んだものです。 そのため、ここでの OUTPUT_POINTS の値自体も、常にこの点数になる事は保証されず、あくまでも目安程度と考えてください。 厳密に保証したい場合は、自由にコードを改造して実装してみてください(これも意外と複雑になります)。

座標値変数の宣言と初期値の設定

続いて、点の動きを追いかけて解曲線を求めていく処理において、点の x, y, z 座標値を格納する変数を宣言しています:

求める解曲線の初期値は、上記変数の初期値として代入しておきます。

解く対象の微分方程式(ローレンツ方程式)の定義

パラメータや変数類の後には、今回数値的に解く対象の方程式であるローレンツ方程式:

\[ \frac{dx}{dt} = - \sigma x + \sigma y \] \[ \frac{dy}{dt} = - x z + \rho x + y \] \[ \frac{dz}{dt} = x y - \beta z \]

の \(\frac{dx}{dt}, \frac{dy}{dt}, \frac{dz}{dt}\)を、それぞれ関数 dxdt, dydt, dzdt として定義しています:

係数が \(\sigma, \rho, \beta\) から S, R, B になっているだけで、内容はそのままですね。

解曲線を求めていく(点の動きを追いかけていく)ループ処理

さて、いよいよ今回のプログラムの本題である、解曲線を求めていく(点の動きを追いかけていく)部分です:

なんだかすごく複雑な雰囲気ですね。 でも、実は処理の流れ自体はシンプルで、 for 文のループ 1 周ごとに、点の座標値 (x,y,z) を時間 DT ぶんずつ動かしていっている だけです。 そのループを何周もぐるぐる回す事で、尺取虫(しゃくとり虫)やパラパラ漫画のように、点が少しずつ動いていくわけです。

そしてループを OUTPUT_INTERVAL 周するごとに、println 関数で x, y, z 座標を一行書き出しているため、 結果のデータには、点の動きの履歴として、解曲線(軌道)が得られるわけです。

点の動きを求める処理のイメージ

さて、ループ1周ぶん、つまり時間 DT ぶんの点の動きは、 x, y, z 座標の時間変化率を与えるローレンツ方程式から、近似的に(ある程度の誤差込みで)求めています。 完璧な精度で求まるとベストなのですが、そのためにはローレンツ方程式を手計算で(解析的に)積分できないといけません。 しかしそれは無理なので、コンピューターで近似的に求めるわけです。

上のコードが複雑になっている理由の大半は、この「 時間 DT ぶんの点の動き 」をなるべく高い精度(少ない誤差で)で求めるために、 4次のルンゲ=クッタ法(RK4)という方法を用いているためです。 多少ラフな精度でよければ、もっとずっと簡単に書けます。 そこで、処理内容の大枠のイメージをつかむために、上のコードよりもずっと単純な方法で「 時間 DT ぶんの点の動き 」を求める事を考えてみましょう。

まず、ある瞬間 \(t\) における点の位置ベクトルを \(\boldsymbol{r_t} = (x_t, y_t, z_t)\) とし、 その瞬間における速度ベクトルを \(\boldsymbol{v_t} = \left( \frac{dx}{dt}, \frac{dy}{dt}, \frac{dz}{dt} \right)_{tの時の値}\) としましょう:

時刻 t における位置ベクトル r と速度ベクトル v

仮に、時間の刻み幅 DT が非常に短い時間なので、その間に点の速度ベクトルはほとんど変化せず、 ほぼ \(\boldsymbol{v_t}\) のままであると近似できるなら、その間は点がほぼ等速直線運動をすると近似できます。 そうすると、時間が DT (数式内では \( \Delta t\) と書く事にします) だけ経過した後の点の位置は:

\[ \boldsymbol{r_{t+\Delta t}} \simeq \boldsymbol{r_t} + \boldsymbol{v_t} \Delta t \]

と近似できます。これをループで繰り返せば、点は少しずつ動いていくはずです。x, y, z 成分ごとにばらしてコードに書くと:

となります(データ出力処理は省いています)。 先ほど見たコードよりも、ずっとシンプルですね。でも目的は同じです。精度が違うだけで、要は DT の間の点の変化を求めているだけ、というわけです。 ちなみに、上の単純な方法は(3次元の)オイラー法といいます。

で、上のオイラー法では、短い時間 DT の間に、速度ベクトル \(\boldsymbol{v}\) がほぼ変化しないと見なしました。 しかし実際にはその間にも微妙に変化するはずで、そのぶんが近似誤差として乗っかってきます。 例えば、電車の速さは1秒の間にはほぼ変化しないでしょうが、1時間の間にはかなり変化しているはずで、後者の場合に速度が一定と近似すると誤差だらけになるでしょう。 つまり速度の変化に対して、時間刻み幅 DT がほどほどに細かくないと、精度よく近似できないわけです。

では、その近似誤差をどうやって減らすかですが、一つは DT をひたすら細かく刻めばいい気がしますね。 でもそうすると、計算の回数がやたらと多くなります。 そして、コンピューターは数値を有限桁で切って扱っているため、基本的には計算を1手行うごとに、わずかな誤差が出ます。 計算の回数があまりにも多いと、そのわずかな誤差が蓄積して、逆に精度が出なくなってしまったりします。 それを置いておいても、計算時間がやたら無駄にかかってしまいます。

という事で「 短い時間 DT の間における \(\boldsymbol{v}\) の変化をもう少し考慮して、DT の刻みが同じでも精度を良くしよう 」という事を考えます。 そのための方法は色々とあり、今回使っている 4次のルンゲ=クッタ法 (RK4) は、その中でも精度と複雑さのバランスがほどよく、結構あちこちで使われるアルゴリズムです。

実際に今回のコード内での該当部分を、データ出力処理を省いてもう一度見てみると:

この通り、色々と中途半端な所を含む地点での \(\boldsymbol{v} = \left( \frac{dx}{dt}, \frac{dy}{dt}, \frac{dz}{dt} \right) \) の値を \(\boldsymbol{k_1}, \boldsymbol{k_2}, \boldsymbol{k_3}, \boldsymbol{k_4}\) として求めて、 それらの重み付き平均のようなものを使う事で、\(\boldsymbol{r} = (x,y,z)\) の変化量を算出している事がわかります。

※ ローレンツ方程式の右辺は時刻 \(t\) を陽には含んでいないため、いまの場合では関数 dxdt, dydt, dzdt の引数に時刻は渡していません。 しかし、時刻を陽に含む微分方程式を扱う場合には、上記コードはもう少し複雑になり、適切に中途半端な時刻値を渡してやる必要があります。 具体的には、k1系には「 t 」、k2系とk3系には「 t + DT/2.0 」、k4系には「 t + DT 」を渡します。

上のように中途半端な点の位置や、謎の重み係数がどこから出てきているのか、また何故それで精度が良くなるのかが気になるかもしれません。 しかし、それを短めでわかりやすく説明する事は少し難しく、恐らくがっつり1回分くらいの記事になってしまいます。 なのでここでは割愛し、またいつかの回に譲る事にします。

その代わりとして、なんとなくの雰囲気だけでも味わいたいという方は、ぜひ併せて以下の記事を読んでみてください:

上の記事では、単純に関数 \(f(x)\) を \(x\) で積分する処理を扱っていますが、精度の異なる複数のアルゴリズム(矩形法、台形法、シンプソン法)を用いて比較しています。 第一歩としては、短い \(\Delta x\) の区間内で \(f(x)\) がほぼ一定と見なし、誤差がどの程度かを評価しています。 そしてその後に、\(\Delta x\) 内での \(f(x)\) の変化を考慮するように改良しながら、精度の向上について評価しています。 最後に登場するシンプソン法では、まさに中途半端な位置での \(f(x)\) の値を用いる事で、精度が大きく向上する様子が見られます。 その様子は、オイラー法から2次&4次のルンゲ=クッタ法への拡張にも通じる所があります。

最後に

今回の内容は以上です。

実は、昔書いたローレンツアトラクタ描画用のGUIツールの記事が、読み返すとちょっと色々と言葉足らずだったので、 今回はVnano版として機能を絞って、その代わりローレンツ方程式の性質とコードの解説を深めに掘り下げようと思って書きました。 …のはずが、今度は逆に長くなり過ぎた気がします。 この種のテーマは掘り下げようと思うとどこまでも掘り下げられるので、なかなかシュッとした感じにまとめるのが難しいですね… すみません。

記事の著者

松井 文宏 - Fumihiro Matsui
[ RINEARN代表、応用情報技術者、博士(理学) ]
VCSSL & Vnano やリニアングラフ、RINPn などを開発しています。 ご意見/ご感想/ご質問などは お問合せページ個人Twitter(雑談歓迎) などへ御気軽にお寄せください。

ライセンス

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


スポンサーリンク



Japanese English
[ 前へ | 目次 | 次へ ]
Vnano版 | 積分値のグラフ描画用データを出力するプログラム

数値的に積分を行い、結果の関数をグラフに描くためのデータを出力するコードです。
3Dグラフを回転アニメーションさせるツール

3Dグラフを、Z軸まわりにゆっくりと回転アニメーションさせるツールです。全角度のグラフを、連番の画像ファイルに保存する事もできます。
連番ファイルから3Dグラフをアニメーション描画するツール

フォルダ内の連番データファイルを読み込み、3Dグラフを高速で連続描画して、アニメーションさせるツールです。グラフを連番の画像ファイルに保存する事もできます。
連番ファイルから2Dグラフをアニメーション描画するツール

フォルダ内の連番データファイルを読み込み、2Dグラフを高速で連続描画して、アニメーションさせるツールです。グラフを連番の画像ファイルに保存する事もできます。
z = f(x,y,t) の形の数式を3Dグラフとしてアニメーション描画するツール

入力欄に z = f(x,y,t) の形の数式を入力すると、それを3次元のグラフにアニメーション描画してくれる簡易ツールです。
y = f(x,t) の形の数式を2Dグラフとしてアニメーション描画するツール

入力欄に y = f(x,t) の形の数式を入力すると、それを2次元のグラフにアニメーション描画してくれる簡易ツールです。
t を媒介変数とする x(t), y(t), z(t) の数式を3Dグラフに描画し、アニメーションもできるツール

入力欄に x(t), y(t), z(t) の形の数式を入力すると、それを t を媒介変数として3次元のグラフに描画し、アニメーションもできる簡易ツールです。
t を媒介変数とする x(t), y(t) の数式を2Dグラフに描画し、アニメーションもできるツール

入力欄に x(t), y(t) の形の数式を入力すると、それを t を媒介変数として2次元のグラフに描画し、アニメーションもできる簡易ツールです。
z = f(x,y) の形の数式を3Dグラフとして描画するツール

入力欄に z = f(x,y) の形の数式を入力すると、それを3次元のグラフに描いてくれる簡易ツールです。
y = f(x) の形の数式を2Dグラフとして描画するツール

入力欄に y = f(x) の形の数式を入力すると、それを2次元のグラフに描いてくれる簡易ツールです。
配列を3Dグラフにアニメーションプロットする(曲面/メッシュグラフ)

座標値配列の内容を、3次元の曲面/メッシュグラフに連続でプロットし、アニメーションさせるサンプルプログラムです。
配列を3Dグラフにアニメーションプロットする(点/線グラフ)

座標値配列の内容を、3次元の点/線グラフに連続でプロットし、アニメーションさせるサンプルプログラムです。
配列を2Dグラフにアニメーションプロットする

座標値配列の内容を、2次元グラフに連続でプロットし、アニメーションさせるサンプルプログラムです。
配列を3Dグラフにプロットする(曲面/メッシュグラフ)

座標値配列の内容を、3次元の曲面/メッシュグラフにプロットするサンプルプログラムです。
ファイルを3Dグラフにプロットする(曲面/メッシュグラフ)

座標値ファイルの内容を、3次元の曲面/メッシュグラフにプロットするサンプルプログラムです。
ユーザーが入力した数式を2Dグラフにプロットする

実行時にユーザーが入力した数式の値を、2次元グラフにプロットするサンプルプログラムです。
配列を2Dグラフにプロットする

座標値配列の内容を、2次元グラフにプロットするサンプルプログラムです。
配列を3Dグラフにプロットする(点/線グラフ)

座標値配列の内容を、3次元の点/線グラフにプロットするサンプルプログラムです。
ファイルを3Dグラフにプロットする(点/線グラフ)

座標値ファイルの内容を、3次元の点/線グラフにプロットするサンプルプログラムです。
ファイルを2Dグラフにプロットする

座標値ファイルの内容を、2次元グラフにプロットするサンプルプログラムです。
この階層の目次
[ 前へ | 目次 | 次へ ]
お知らせ

Vnanoがオープンベータ版に移行、VCSSLの実行環境で標準で実行可能に
2021年04月07日 - ソフト内組み込み用スクリプトエンジン&言語「 Vnano 」がオープンベータ版に移行し、併せて、VCSSLの実行環境でもVnanoコードの実行が可能になりました。詳細をお知らせします。

リニアングラフの最新版をリリース、2D版でも描画エンジンの直接操作が可能に
2021年04月03日 - リニアングラフ2D/3Dの最新版をリリースしました。それぞれのアップデート内容をお知らせします。今回から、2D版でもJava言語APIによる描画エンジンの直接操作が可能になりました。

RINPn のオープンベータ版をリリース! 詳細な公式ガイドも同梱&公開
2021年03月08日 - 2019年より開発進行中のプログラム関数電卓「 RINPn(りんぷん)」が、正式リリースに向けた最終準備段階として、オープンベータ版へと移行しました。その詳細をお知らせします。

新着
Vnano版 | ローレンツ方程式を数値的に解くプログラム

ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。
2021年02月12日
Vnano版 | 積分値のグラフ描画用データを出力するプログラム

数値的に積分を行い、結果の関数をグラフに描くためのデータを出力するコードです。
2020年12月20日
Vnano版 | 積分値を求めるプログラム (数値積分)

矩形法/台形法/シンプソン法を用いて、積分の値を数値的に求めるコードです。
2020年12月20日
3Dグラフを回転アニメーションさせるツール

3Dグラフを、Z軸まわりにゆっくりと回転アニメーションさせるツールです。全角度のグラフを、連番の画像ファイルに保存する事もできます。
2019年10月09日
[公式ガイドサンプル] ユーザーのGUI操作に対して処理を行う

「VCSSL GUI開発ガイド」内のサンプルコードです。ユーザーがGUIを操作した際に行う処理を実装します。
2019年07月28日
開発元Twitterアカウント