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

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

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

スポンサーリンク

---

使用方法

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

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

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

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

プログラムの起動

Windows をご使用の場合

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

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

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

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

Linux 等をご使用の場合

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

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

起動後

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

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

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

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

続いて、ここで表示/保存したデータからグラフを描いてみましょう。すぐ下の節で説明します。

グラフの描画方法

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

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

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

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

このグラフが、プログラム内で指定されたパラメータ・初期値の条件下で、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項目に、このファイル名をそのまま入力します。 すると黒い画面が表示され、そこにグラフ用のデータが出力されます：

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

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

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

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

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

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

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

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

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

---

スポンサーリンク

---

題材解説

ローレンツ方程式

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

ローレンツ方程式は、アメリカの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 つの点の動きを、それぞれ赤色と緑色の点として、重ねてアニメーションさせているものです。

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

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

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

予測の難しさ

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

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

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

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

コード解説

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

コード全体

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

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

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

大体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の時の値}\) としましょう：

仮に、時間の刻み幅 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)\) の変化量を算出している事がわかります。

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

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

• 積分値を求めるプログラム (数値積分)

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

最後に

今回の内容は以上です。

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

ライセンス

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