ExODEPI - デカルト言語 - 並列ロジック プログラミング言語 Wiki - デカルト言語 - 並列ロジック プログラミング言語

←　前のページに戻る

2.1 π(パイ)

これは、微分方程式とはちょっと違うのですが、ODE機能の簡単な使い方を示す最初の例題として選びました。

ここでは数値積分を行わせます。微分方程式を数値計算して解くためには、数値積分が必要だからです。

さて、以下のような積分を計算しましょう。

　　　　　1
  y = ∫ 4/(1+x^2) dx
　　　　　0

これは4/(1+x^2)を0から1の範囲で積分します。

ここでは数学的な説明はしませんが、答は”π　”(3.14159...)になります。

このような積分の計算をまずデカルト言語のODE機能で行わせてみましょう。

2.1.1 pai.dec ソース

/* π */

? ::sys <PrintResultOff>
  ::sys <ODE (#x 0 1 0.001)
        <letf #y = 4 / (1 + #x * #x)>
        ::sys <integral #y2 #y>
        ::sys <ODEprint 100 <%_"1.1f" #x>  "," <%_"1.17f" #y2>>
  >
  ;

2.1.2 ソース解説

上記のソースについて解説します。

まず、最初のPrintResultOffというのは、余分な結果出力を抑止するためのものです。

次の行のODE述語というのが、ODE機能の中心の述語です。

この述語のマニュアルを次に示します。

::sys <ODE [VAR] (VAR INIT END STEP) PRED-LIST>

        微分方程式(ODE)を解く述語。
        初期値INITからSTEP単位でENDまでの値が変数VARに設定
        されて繰り返しPRED-LISTの術語を実行する。

例では、ODE述語の第1引数が(#x 0 1 0.001)となっているので、#xという変数に0から1までの間を0.001刻みで繰り返すことを意味します。

これだけを見ると単なる繰り返しを行うだけに見えますが、実はそれだけではありません。 この述語の第2引数から指定される述語には他の場所では使えないODE関連の述語を指定することができるのです。 表からは見えないODE関連の裏の処理をサポートしている述語でもあります。

ODE述語の中に上の例では３つの述語があります。 これらの第2引数から含まれている述語を順に実行していきます。

まず最初の<letf #y = 4 / (1 + #x * #x)>は、単純な代入式です。y = 4 / (1+x^2)を計算しているのですね。

次の::sys <integral #y2 #y>が、今回の処理の肝となる処理です。 ODE述語の中でintegral述語を使うと、第2引数に指定した値の積分値を第1引数の変数に設定します。 積分の計算には、Simpson法を使っています。

integral述語は数値積分を行うためにODE述語の内部処理を利用しているため、ODE述語の外では使うことができない述語です。

３つ目の::sys <ODEprint 100 <%_"1.1f" #x> "," <%_"1.17f" #y2>>も同様にODE述語の中でのみ使える述語です。 第1引数の100というのは、100回実行されるごとに第2引数以降の値を表示する述語です。 <%_"1.1f" #x>は、変数#xの値をC言語の"%1.1f"と同じフォーマットで表示します。 ここでは、0から1までの#xと積分値#y2を、100個おきに表示していきます。

なぜ、ODEprintを使うかというと、次に示すような理由です。

数値積分は細かな刻みで繰り返さないと精度が得られません。 しかし、だからといって結果をすべて表示するとあまりに多数の結果が得られてしまいます。 そのような場合に結果を間引いて適切な刻み幅で表示したいようなときに使います。

2.1.2 実行結果

では、1.1.1で示したソースをpai.decという名前で保存して実行してみましょう。

$ descartes pai.dec
 0.0 , 0.00000000000000000
 0.1 , 0.39867460996469886
 0.2 , 0.78958223939961058
 0.3 , 1.16582717791157154
 0.4 , 1.52202550844955854
 0.5 , 1.85459043600330642
 0.6 , 2.16167800108239652
 0.7 , 2.44290385755687313
 0.8 , 2.69896376889423192
 0.9 , 2.93126040714603490
 1.0 , 3.14159265358979328

","(カンマ)で区切られた２列の数字列が表示されます。 順に0からカンマの左に表示された値までの積分値が、カンマの右に表示されます。

最終的に最後に表示された「 1.0 , 3.14159265358979328 」が0から1まで積分した値です。

完全なπ（円周率 3.14159265358979323846...）とはちょっと違いますが、結構よい精度で計算できていますね。

2.1.3 グラフ

先に示した結果出力は、数値を","(カンマ)で区切った形式で出力していました。 このような、カンマで区切られた数値列は、EXCELなどの表計算のプログラムではcsv(Comma Separated Values：カンマ区切り)形式と呼ばれ、そのまま読み込んで処理することができます。

以下のように、pai.csvという名前で結果を保存してみましょう。

$ descartes pai.dec > pai.csv

このように保存されたファイルは表計算ソフトで使えます。 Windowsならばpai.csvをダブルクリックしただけでEXCELで開くことができるはずです。

試しにEXCELで読み込んでから、グラフ化したものを次に示します。

(グラフ化する方法は、お持ちの表計算プログラムの方法に従ってください。)