GNU Octave 2.1.x 日本語マニュアル/数値積分

22　積分　Quadrature

22.1 一変数関数

[v, ier, nfun, err] = quad (f, a, b, tol, sing)

[Loadable Function]

QUADPACKを用いて一変数関数の非線形数値積分を行います。 最初の引数は関数の名前で積分を計算する関数のハンドルあるいはインライン関数です。関数は以下のフォームをもたなければなりません。

$y=f(x)$

ここで　x と y はスカラーです。 二番目と三番目の引数は、積分範囲です。 片方あるいは両方の値は無限大になり得ます。 オプションの引数tolはベクトルで結果の希望する精度を指定します。 ベクトルの最初のエレメントは希望する絶対精度で、二番目のエレメントは希望する相対精度です。 相対精度テストだけ行うとき絶対精度はゼロにおきます。 絶対精度でテストだけ行うとき相対精度はゼロにおきます。 オプショナル引数singは、積分が特異の時の数値ベクトルです。 積分の結果はｖに返されます、ierは整数のエラーコードです（0は正常に積分が行われたことを示します）。 nfunは、何回関数が評価されたかを示します。errは、解の誤差の推定値を含みます。

quad_options (opt, val)

[Loadable Function]

二つの引数で呼び出されると、この関数は、関数のquadのオプションパラメータを設定できます。 一引数が与えらるれと、関数quad_optionsは、対応するオプションの値を返します。 引数を指定しない場合、利用可能なすべてのオプションとその現在の値の名前が表示されます。

"absolute tolerance":"絶対許容誤差"

絶対許容誤差は、純粋な相対的なエラー・テストのばあいはゼロとなります。

"relative tolerance":"相対的許容値"

非負の相対許容誤差。

絶対許容誤差(Absolute tolerance)がゼロであれば、相対的許容値誤差(Relative　tolerance)は,

$max(50*EPS,0.5e-28)$

より大きいか、または等しくなければなりません。

これは、quadを使った関数の積分例です。

$y=x\sin(1/x){\sqrt {|{1-x}|}}$

from x = 0 to x = 3.

これは、かなり難しい積分です（難しい理由を確かめるために積分範囲をプロットしてください）. 最初のステップは関数を定義します：

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

演算子‘dot’フォームの使用に留意してください。quadの呼び出しで必須ではありませんが、プロットのための点の集合を生成するのがはるかに簡単になります（引数がベクトルで、結果がベクトルを持つ関数として呼び出すことを可能にします）。 単にquadを呼び出します:

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

難しい積分にもかかわらず、quadがierにゼロ以外の値を返したら結果は十分に正確です。 （理由を理解するために、下限を0.1,　次に 0.01,　さらに 0.001, etc.と移動したときに何がおこるか試験してください）

22.2 Orthogonal Collocation

[r, amat, bmat, q] = colloc (n, "left", "right")

[Loadable Function]

直交配列法の微分と積分の重み行列を、J. Villadsen and M. L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation 、で与えられたサブルーチンを使って計算します。 二階微分方程式をcollocの重み行列を生成して解く例題を示します。

$u'-{\alpha }{u''}=0$

境界条件

$u(0)=0$ and $u(1)=1$

まず、n個の点（区間の端点を含む）のための重み行列を生成し、右辺の境界条件を（$\alpha$の特定の値の場合）に組み込むことができます。

n = 7;
alpha = 0.1;
[r, a, b] = colloc (n-2, "left", "right");
at = a(2:n-1,2:n-1);
bt = b(2:n-1,2:n-1);
rhs = alpha * b(2:n-1,n) - a(2:n-1,n);

根rにおける解は、

u = [ 0; (at - alpha * bt) \ rhs; 1]
) [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]