求積

Octave では Quadpack を使って、1変数非線型関数の積分を計算することができます。積分計算の関数 quad の書式は次のようになります。

[V, IER, NFUN, ERR] = quad (F, A, B, TOL, SING)

最初の引数の F は積分をする関数名です。積分する関数は次の型をしている必要があります。

y = f (x)

この式の変数の y と x はスカラーです。

2 番目と 3 番目の引数 A と B は積分区間です。どちらも無限大を指定することができます。

省略可能な引数の TOL はベクトルのデータで、積分結果の精度を指定することができます。ベクトルの一番目の要素で、絶対誤差の許容範囲を指定し、二番目の要素で相対誤差の許容範囲を指定します。相対誤差のみで積分計算を終了させる場合、絶対誤差範囲を 0 にしておきます。絶対誤差のみを採用する場合、相対誤差の方を 0 にします。

省略可能な引数 SING は数値ベクトルで、その点で被積分関数の値が単一であることが分かっている場合に指定します。

積分の結果は戻り値 V に渡され、積分時のエラーコードは IER に渡されます。0 が積分が成功したときのコードです。NFN は 被積分関数の評価が何回必要だったかを示し、ERR には解の推定誤差が入ります。

quad 関数のオプションを取り扱うのが次の quad_options 関数です。

quad_options (OPT, VAL)

OPT と VAL の二つの引数が指定されると、OPT に指定された quad のオプションに VALの値がセットされます。引数が OPT だけの場合 OPT のオプションの値が表示されます。引数なしで quad_options 関数を呼び出すと、全てのオプションの名前と値が一覧表示されます。

quad を利用して次の関数の数値積分を行ってみましょう。

F(X) = X * sin (1/X) * sqrt (abs (1 - X))

ここで積分区間は X = 0 から X = 3 までとします。

これは割りに難しい積分です。(gnuplot の plot コマンドで次のようにして確認できます。)

gnuplot> set grid
gnuplot> plot [-1:4][:] x * sin (1/x) * sqrt (abs(1 - x))

quad による計算の最初のステップは関数を定義することです。

     function y = f (x)
       y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
     endfunction

関数の定義にドットつき演算子を用いていますが、こうすることによって Y = f (X) のような演算が引数や戻り値がベクトルであっても計算することができます。関数の定義が済めば後は quad を呼ぶだけで積分を計算することができます。ただし、関数の定義はドット演算子ではないふつうの演算子を使っても行うことができます。しかしその場合は引数 x はスカラーに限定されます。

     [v, ier, nfun, err] = quad ("f", 0, 3)
          => 1.9819
          => 1
          => 5061
          => 1.1522e-07

IER の値が 0 ではありませんが err の数値から見ると充分に精確な結果と言えます。

この項の説明は丸々マニュアルの翻訳です。