9.2 変数変換型積分公式, 適応自動積分
9.2.1 台形則の誤差
台形則では分点数を倍 (刻み幅 h を 1/2) にすると誤差はおおよそ 1/4 になった. これは, 被積分関数が解析的な場合台形則の誤差が次のように表されることによる (オイラー・マクローリン展開).
\[
I_m – I = c_2 h^2(f'(b) – f'(a)) + c_4 h^4(f”'(b) – f”'(a)) + \dots + c_2p h^{2p}(f^{(2p-1)}(b) – f^{(2p-1)}(a)) + O(h^{2p+2})
\]
ここで \(c_i\) は \(h\) や \(f\) に無関係な定数である.
\(h^4\) 以上の項が十分小さいと仮定すると誤差はほぼ \(h^2\) に比例すると見なすことができ, \(h\) が 1/2 になると誤差は 1/4 になる. しかし, 微分不可能な特異点があったりするとこれが成り立たなくなることは「9.1.3 数値実験」で見た.
数値実験 (7) 周期関数
特別な場合の例として次の関数の [0, \(\pi\)/2] での積分を求める.
\[
I = \int_0^{\pi/2} \frac{1}{\sqrt{1 – sin^2(x) / 2}} dx = 1.85407467730137
\]
図からもわかるようにこれは周期関数である. 種々の積分公式で分点数 N (横軸) を変えて計算し, その時の相対誤差 (縦軸) を対数グラフにすると次のようになった.
意外にも台形則が最も速く収束している. 台形則の結果を表にしてみると, 収束率 = 1/4 と比較すると桁違いに小さく (収束が速く) なっていることがわかる. こうなると図のようにロンバーグ積分はかえって遅くなる.
| 分点数 | 積分値(計算値) | 積分値(真値) | 誤差 | 収束率 |
|---|---|---|---|---|
| 2 | 1.85495913108563 | 1.85407467730137 | 8.8445378e-4 | |
| 4 | 1.85407522776731 | 1.85407467730137 | 5.5046594e-7 | 1/1607 |
| 8 | 1.85407467730167 | 1.85407467730137 | 2.9465319e-13 | 1/1.868e+6 |
| 16 | 1.85407467730137 | 1.85407467730137 | 0 |
数値実験 (8) 無限積分
もうひとつ, 無限積分の例を見てみる.
\[
\int_0^{\infty} e^{-x^2/2} dx = 1.25331413731550
\]
関数形は次のとおりである. 関数値は x = 8 で \(10^{-14}\) のオーダー, x = 10 で \(10^{-22}\) のオーダーと急速に 0 に近づいていく関数になっている.
ここで, 積分範囲を \([0, \infty]\) の代わりに \([0, 10]\) として計算してみると次のようになった.
ここでも台形則が一番速く収束し, N = 15 で最高精度に達した.
減衰が速い関数では端ですぐに微分係数が 0 に近づくため, 上式で表される誤差も速く 0 に近づくと考えられる. このような無限積分に対しても台形則が他の公式よりも精度が良いことが知られている.
9.2.2 変数変換型積分公式
台形則は周期性のある関数の積分や減衰の速い無限積分に対して有効であることがわかったが, 実際に台形則がそのまま使われることはない. しかし, 台形則のこの特長は以下に示す変数変換による積分公式に応用されている.
一般の積分を「9.2.1 台形則の誤差」の 2 つのケースのように台形則が速くなるような積分になるように変数変換し, それに台形則を適用すれば性能の良い積分公式となることが期待できる. このアイデアに基づく公式が変数変換型積分公式である.
積分変数を変換する公式は次のとおりである.
\[
\begin{align}
& \int_a^b f(x)dx = \int_c^d g(t)dt \\
& ただし, \space x = \phi(t), \space g(t) = f(\phi(t))\phi'(t) \\
\end{align}
\]
ここで \(\phi(t)\) は区間 [c, d] で単調で, かつ \(a = \phi(c), \space b = \phi(d)\) である. これにより積分値が不変のまま被積分関数を変換することができる. 適切な \(\phi(t)\) を見つけることができれば効果的な公式を作ることができる.
9.2.2.1 IMT公式
\(\phi(t)\) として次式を使うのが IMT 公式である.
\[
\begin{align}
& \phi(t) = \frac{1}{Q} \int_0^t \psi(r) dr, \space Q = \int_0^1 \psi(r) dr \\
& ただし, \space \psi(r) = exp(\frac{-1}{r} – \frac{1}{1 – r}) \\
\end{align}
\]
この変換を区間 [0, 1] の積分に適用すると, 変換後の関数 \(g(t)\) は \(g^{(j)}(0) = g^{(j)}(1) = 0 \space (j = 0, 1, \dots)\) となる. すなわち, 積分区間の両端で微分係数が 0 になっていて, 台形則の誤差が 0 に近くなる条件を満たしている.
この変換を施し台形則を適用すると, 分点数を N として次の計算式が得られる.
\[
\begin{align}
& S_N = \frac{1}{N} \sum_{i=1}^N w_i^{(N)} f(x_i^{(N)}) \\
& 重み: w_i^{(N)} = \phi'(\frac{i}{N}) = \frac{1}{Q}\psi(\frac{i}{N}) \\
& 分点: x_i^{(N)} = \phi(\frac{i}{N}) = \frac{1}{Q}\int_0^{\frac{i}{N}} \psi(r) dr \\
\end{align}
\]
IMT 公式については文献[3], [4]に詳しい.
9.2.2.2 二重指数型積分公式 (DE 公式)
対象とする積分を減衰の速い無限積分に変換して台形則の誤差が 0 に近くなる条件を満たすことにより効率よく積分を計算するのが二重指数型積分公式 (DE 公式) である. 変換関数 \(\phi(t)\) は積分区間などによっていくつかの種類が提案されている.
DE 公式については文献[2], [4], [5]に詳しい. 変数変換型積分公式としては IMT 公式がはじめに提案され, その後 DE 公式が提案された. 現在では主に DE 公式が使われている.
(1) 有限区間 [-1, 1] の積分
\(\phi(t)\) として次を使う.
\[
\begin{align}
& S_h^{(N)} = h \sum_{i=N_-}^{N_+} w_i^{(N)} f(x_i^{(N)}) \\
& \phi(t) = tanh(\frac{\pi}{2}sinh(t)) \\
& \phi'(t) = \frac{\frac{\pi}{2}cosh(t)}{cosh^2(\frac{\pi}{2}sinh(t))} \\
\end{align}
\]
きざみを \(h\), 分点数を \(N (= N_+ + N_- + 1)\) とすると, 計算式は次のようになる.
\[
\begin{align}
& 重み: w_i^{(N)} = h \phi'(ih) \\
& 分点: x_i^{(N)} = \phi(ih) \\
\end{align}
\]
無限積分の計算は下側を \(-N_-\), 上側を \(N_+\) で打ち切る. これらの値は許容誤差限界をもとに適切な値を設定する.
端点に特異点がある関数の例で取り上げた \(I = \int_0^1 \sqrt{1 – x^2} dx = \pi/4\) の関数のグラフを再掲する. この関数は下の数値実験 (2) でも使用する.
これを上の変換関数 \(\phi(t)\) を使って無限積分に変換すると次のようになる.
上側も下側も減衰が速いことがわかる. きざみ幅を 0.1 としてやると, 上側・下側それぞれ 32 点で \(10^{-15}\) 以下のオーダーまで減衰した. これは台形則で精度よく積分することができる.
(2) 半無限区間 [0, \(\infty\)] の積分
半無限区間の積分において被積分関数がゆるやかに減衰するような場合についても有限区間の場合と同様の考え方で DE 公式を導くことができる. 文献[2]に示されているプログラム DEHINT では被積分関数により 3 種類の公式を使い分けるようになっている.
(a) 通常の関数の場合
\[
\begin{align}
& \phi(t) = exp(2sinh(t)) \\
& \phi'(t) = 2cosh(t)exp(2sinh(t)) \\
\end{align}
\]
(b) 被積分関数が指数関数的に減衰する場合 (\(f(x) = f_1(x)exp(-x)\) の形をしている場合)
\[
\begin{align}
& \phi(t) = exp(t – e^{-t}) \\
& \phi'(t) = (1 + e^{-t})exp(t – e^{-t}) \\
\end{align}
\]
(c) 被積分関数の減衰がさらに急激な場合 (\(f(x) = f_2(x)exp(-x^2)\) の形をしている場合)
\[
\begin{align}
& \phi(t) = exp(\frac{t}{2} – e^{-t}) \\
& \phi'(t) = (\frac{1}{2} + e^{-t})exp(\frac{t}{2} – e^{-t}) \\
\end{align}
\]
(3) 全無限区間 [\(-\infty\), \(+\infty\)] の積分
被積分関数がゆるやかに減衰するような全無限区間の積分についても同様に DE 公式を導くことができる.
\[
\begin{align}
& \phi(t) = sinh(\frac{\pi}{2}sinh(t)) \\
& \phi'(t) = \frac{\pi}{2}cosh(t)cosh(\frac{\pi}{2}sinh(t)) \\
\end{align}
\]
(4) フーリエ型積分
次の半無限区間の積分 (フーリエ型積分) を求める.
\[
I_c = \int_a^{\infty} f(x) cos(\omega x) dx \space または \space I_s = \int_a^{\infty} f(x) sin(\omega x) dx
\]
次の変換関数を使用して無限区間 [\(-\infty\), \(+\infty\)] の積分に変換する.
\[
\begin{align}
& I_c の場合: x = M \phi(t – \frac{\pi}{2M})/\omega
& I_s の場合: x = M \phi(t)/\omega
\end{align}
\]
\(\phi(t)\) は \(t \to -\infty\) のときに二重指数関数的に \(\phi'(t) \to 0, \space t \to +\infty\) のときに二重指数関数的に \(\phi(t) \to t\) となるような関数である. きざみ幅 \(h\) と定数 \(M\) は, \(Mh = \pi\) が成り立つように選ばれる. この変換により, x が大きくなるときに公式の分点が \(cos(\omega x)\) あるいは \(sin(\omega x)\) のゼロ点に二重指数関数的に近づくようになり, 大きな x については関数値を計算しなくてもよくなる.
φ(t)はいくつか提案されているが, 例えば次式がある.
\[
\begin{align}
& \phi(t) = \frac{t}{1 – exp(-ksinh(t))}, \space k = 6 \space (または \space 2\pi) \space (文献[4], [5], [7]) \\
& または \\
& \phi(t) = \frac{t}{1 – exp(-2t -\alpha(1 – e^{-t}) – \beta(e^t – 1))}, \space \beta = \frac{1}{4}, \space \alpha = \frac{\beta}{\sqrt{1 + M log(1 + M)/4\pi}} \space (文献[9]) \\
\end{align}
\]
下の数値実験 (3) で使用する \(I = \int_0^{\infty} log(x)sin(x)dx = -\gamma\) の関数のグラフを示す.
文献[7]の公式 (\(K = 2\pi\)) を使い, \(Tol = 10^{-5}\) としたときの分点を赤丸で示した. x が大きくなると分点が関数のゼロ点に一致していくのがわかる.
数値実験 (9) 性質のよい関数
「9.1.3.1 (1) 性質の良い関数」で実験した関数を用いて変数変換型公式の性能を見る.
次の積分を計算する。
\[
I = \int_0^1 e^x cos(x) dx = 1.378024613547364
\]
DE 公式ではプログラムの作りによって関数評価回数が変わってくるが, ここではきざみみ幅を変化させてそのときの実際の精度にみあったところで無限積分を打ち切るようにした (自動積分ではない).
比較のため台形則, シンプソン則, ガウス公式もプロットした.
IMT 公式はほぼシンプソン則なみの性能に見える. DE 公式は IMT 公式よりは速いがガウス公式ほどではない.
数値実験 (10) 端点に特異点を持つ関数
「9.1.3.1 (3) 端点に特異点を持つ関数」で実験した端点に特異点を持つ関数の場合について見る.
\[
I = \int_0^1 \sqrt{1 – x^2} dx = \frac{\pi}{4}
\]
x = 1 に特異点がある.
前に示したように補間型公式ではどれも十分な結果が得られなかったが, IMT 公式と DE 公式は特異点の影響を受けずに性質のよい関数のときとほぼ同じように計算できた. このように, 端点に特異点がある (端点で微分係数が定義できない) 場合にも適用できることは IMT 公式と DE 公式の特長のひとつである.
数値実験 (11) フーリエ型積分
次のフーリエ型積分を計算する.
\[
I = \int_0^{\infty} log(x)sin(x) dx = -\gamma
\]
ここでは文献[7]の公式を使用した (\(K = 2\pi\)). 要求精度 Tol を適当に変化させ, \(M = -\pi log(Tol)/2\) と固定して計算した. ちなみに, XLPack に収録されている Deoint では推定誤差をみて M の値を調節するようになっている.
9.2.3 誤差の推定
「9.1.3 数値実験」で見たように関数評価回数が同程度であれば種々の公式の中でガウス公式が最も精度が良かった. 単に計算するのであればこれで良いが, 実務的には積分値と同時に推定誤差も知りたい. また, 後で示す自動積分のベースとして使う場合には誤差推定ができなければならない. ガウス公式の弱点は誤差の推定がしにくいことである.
9.2.3.1 台形則の誤差推定
台形則の計算結果の誤差を推定するためには, 例えば, N 分点で計算したものと 2N 分点で計算したものを比較するなどすればよい. 次の例 (滑らかな関数) で台形則の計算例を見る.
\[
I = \int_0^1 e^x cos x dx = 1.37802461354736
\]
| 分点数 | 積分値(計算値) | 積分値(真値) | 誤差 | 誤差推定値 |
|---|---|---|---|---|
| 2 | 1.34061800327106 | 1.37802461354736 | 0.037406610 | 0.035423678 |
| 4 | 1.36858238253106 | 1.37802461354736 | 0.009442231 | 0.009321460 |
| 8 | 1.37565843490021 | 1.37802461354736 | 0.002366179 | 0.002358684 |
| 16 | 1.37743271822098 | 1.37802461354736 | 0.000591895 | 0.000591428 |
| 32 | 1.37787661780930 | 1.37802461354736 | 0.000147996 | 0.000147967 |
\[
推定誤差 = \frac{4I – I_{prev}}{3}
\] ここでは分点数を倍々にするような順序 (等間隔な分割幅それぞれを半分にしていく) で計算を行った. \(I\) は積分の計算値, \(I_{prev}\) は分点数が半分の時の (つまり 1 つ前の) 計算値である. この場合, 分点数を倍にした時にもう一回最初から全部の分点で計算する必要はなく, 半分はすでに計算済の値がそのまま使える.
9.2.3.2 誤差推定機能付きのガウス公式
台形則とは異なり, ガウス公式では分点が等間隔でないため, 分点数を倍にするともう一回最初から全部の分点で計算する必要があるのが弱点である.
これを解決するために, N 点ガウス公式の計算値を生かしたまま 2N + 1 点を使った積分公式がガウス・クロンロッド公式である.
N 点ガウス公式による積分を \(G_N\) で表す (積分区間を [-1, 1] とする).
\[
G_N = \sum_{i=1}^N w_i f(x_i)
\]
ガウス・クロンロッド公式では同じ N 個の分点 \(x_i\) と新たに計算が必要な N + 1 個の分点 \(y_i\) を使って求める積分 \(K_{2N+1}\) を次式で表す.
\[
K_{2N+1} = \sum_{i=1}^N a_i f(x_i) + \sum_{j=1}^{N+1} b_j f(y_j)
\]
N 個の分点 \(x_i\) における係数 \(a_i\) はガウス公式の係数 \(w_i\) とは異なる値になる.
N 点ガウス公式 \(G_N\) は 2N – 1 次多項式まで正しい解を与える (すなわち 2N – 1 次である) が, ガウス・クロンロッド公式 \(K_{2N+1}\) は 3N + 1 次 (N が偶数のときは 3N + 2 次) となるように分点と重みを選ぶ.
\(x_i\) は N 次のルジャンドル多項式 \(P_N(x)\) のゼロ点であるが, 残りの分点 \(y_i\) が n + 1 次の多項式 \(E_{N+1}(x)\) のゼロ点になるとすると, \(E_{N+1}(x)\) は次の直交関係を満たす.
\[
\int P_N(x) E_{N+1}(x) x^k dx = 0 \space (k = 0, \dots, N)
\]
上の条件を満たす多項式 \(E_{N+1}(x)\) のゼロ点は連立方程式を解いて求めることができ, 区間 (-1, 1) に存在し, \(P_N(x)\) のゼロ点と交互に現れることがわかっている. クロンロッドはこの連立方程式を解いて求めたが, この方程式は非常に悪条件で解きにくいため, 後にいくつかの計算法が提案されている.
N = 2 ～ 10 のときの分点および重みの値を「9.2.5 積分公式の係数表」に掲載した. ここでは, 文献[8]に従い, ヤコビ・クロンロッド行列とよばれる対称三重対角行列を作りその固有値を求める方法により計算した.
ガウス・クロンロッド公式では, 2N + 1 回の関数評価で 2N + 1 次の \(G_N\) と \(3N + 1\) 次の \(K_{2N+1}\) の両方を求めることができ, この 2 つを比べることにより誤差を推定することができる. 推定誤差は経験的に次式で表される.
\[
推定誤差 = (200*|G_N – K_{2N+1}|)^{1.5}
\]
この式による推定誤差は安全側に求められるとされている.
ガウス・クロンロッド公式 \(K_{2N+1}\) は 3N + 1 次であるが, 同じ関数評価回数で計算できるガウス公式 \(G_{2N+1}\) は 4N + 1 次である. ガウス・クロンロッド公式は精度的にはガウス公式に劣るが, 誤差の推定ができる利点の方が大きい場合に有効である.
数値実験 (12) ガウス・クロンロッド公式
ガウス・クロンロッド公式を用いて次の積分を求める.
\[
I = \int_0^1 e^x cos(x) dx = 1.37802461354736
\]
これは前にも例題に使用した性質のよい関数であり, 計算結果は次のようになった.
横軸に関数評価回数, 縦軸に誤差をプロットした. 緑色のデータは前項の式で求めた推定誤差である. 横軸は分点数 (ガウス公式では N, ガウス・クロンロッド公式では 2N + 1 の値) に等しい.
推定誤差はかなり安全側に算出されているのがわかる.
数値実験 (13) ガウス・クロンロッド公式 (特異性のある場合)
上の例題では収束が速くすぐに丸め誤差の領域に入っていてうまく比較できないのでもう少し難しい次の積分を計算する.
\[
I = \int_0^4 1/(1 + x)dx = ln(1 + 4) = ln(5)
\]
これも前に例題に使用した関数で, 積分区間外であるが近くに特異点を持つ.
この例ではガウス・クロンロッド公式を使ってもそれほど大きなペナルティを払わずに誤差を推定できている.
9.2.4 自動積分プログラム
自動積分プログラムとは, 積分区間と要求精度を与えると誤差推定値が要求精度を満たすまできざみ幅を順次小さくするなどして積分を繰り返すプログラムである.
適応自動積分プログラムはこれを一歩進めて, 積分区間を部分区間に分割し, 誤差推定値が要求された精度を満たすように部分区間の幅を自動的に調節するものである. 被積分関数がなめらかにゆっくり変化するような区間では大きな区間幅を用い, 誤差推定値が大きい部分では細かい区間に分割してきざみ幅を小さくして計算することにより, 最小の計算量で要求精度の結果を得ようとするものである.
適応自動積分としては, 誤差推定を行うためと計算済の関数値を再利用するために, 分点が等間隔なシンプソン則やニュートン・コーツ公式が使われてきた. しかし, 現在では精度を追求して, 誤差推定機能を持つガウス・クロンロッド公式のような不等間隔分点ではあるが高精度の公式を使ったものが使われている.
9.2.4.1 適応自動積分プログラム QAG
適応自動積分ルーチンの例として QUADPACK (文献[6]) の QAG について簡単に説明する. QAG は要求精度を満たすまで自動的に積分区間を部分区間に分割しながらガウス・クロンロッド公式の固定分点プログラム (QK15, …, QK61) を呼び出すことにより自動積分を行う. XLPack には Qag として収録されており, 通常の数値積分プログラムとしては最も頻繁に使用される.
このプログラムの動作の概要は次のとおりである.
与えられた収束条件を「推定誤差 \(\le\) Tol」とする.
- まず, 全区間をガウス・クロンロッド公式の固定分点プログラムにより積分を行い積分値 S と推定誤差 E を求める.
- E \(>\) Tol である限り以下を繰り返す.
- 推定誤差が最も大きい部分区間を半分に分割してそれぞれを固定分点プログラムにより積分し, 積分値 \(S_i\) と推定誤差 \(E_i\) を求める(i は部分区間ごとの値であることを表す).
- 新たに構成されたものを含めたすべての部分区間における積分値と推定誤差のそれぞれの総和 S と E を更新する.
- E \(\le\) Tol になったらその時の積分値の総和 S を最終結果として返す.
QAG ではベースの固定分点プログラムとして QK15, …, QK61 のうちユーザーが指定したものを使用し, これを自動的に変更することはしない. また, 区間を分割する際に計算済の関数値を再利用することもしない.
数値実験 (14) Qag (Qk15)
例として, 区間 [0,1] において次の関数 (文献[1]参照) を積分した時の QAG の挙動を下図に示す.
\[
f(x) = \frac{1}{(x – 0.3)^2 + 0.01} + \frac{1}{(x – 0.9)^2 + 0.04} – 6
\]
ここでは, XLPack に収録されている Qag を使用し, 目的誤差を \(10^{-5}\) として固定分点プログラムは Qk15 を指定して計算した. 区間は 5 つに分割された (0 ～ 0.25, 0.25 ～ 0.375, 0.375 ～ 0.5, 0.5 ～ 0.75, 0.75 ～ 1). トータルの関数評価回数は 135 回であった. 最も急峻な変化のところで分割されているのがわかる.
数値実験 (15) Qag (Qk61)
同じ条件で Qk61 を指定すると次のようになった.
2 分割で計算しており (0 ～ 0.5, 0.5 ～ 1), トータルの関数評価回数は 183 回になった. 変化が激しい場合には QK15 を選択し分割数を多くし, そうでない場合には QK61 のように高精度のものを選び分割数を少なくするとよさそうである.
9.2.5 積分公式の係数表
(1) ガウス・クロンロッド公式
\[
K_{2N+1} = \sum_{i=1}^N a_i f(x_i) + \sum_{j=1}^{n+1} b_j f(y_j)
\]
ガウス・クロンロッド公式における分点 \(x_i, y_j および重み a_i, b_j\) を以下に示す(\(x_i\) はガウス公式の \(x_i\) と同じ). 計算は拡張倍精度 (80 ビット浮動小数) により行い, 結果の 16 桁を表示した.
N は対応するガウス公式の N の値を示す. ガウス・クロンロッド公式としての分点数は 2N + 1 となる (例. N = 7 は 15 分点).
N = 2
xi: -5.773502691896258E-0001 5.773502691896258E-0001
ai: 4.909090909090909E-0001 4.909090909090909E-0001
yi: -9.258200997725515E-0001 0.000000000000000E+0000 9.258200997725515E-0001
bi: 1.979797979797980E-0001 6.222222222222222E-0001 1.979797979797980E-0001
N = 3
xi: -7.745966692414834E-0001 0.000000000000000E+0000 7.745966692414834E-0001
ai: 2.684880898683334E-0001 4.509165386584741E-0001 2.684880898683334E-0001
yi: -9.604912687080203E-0001 -4.342437493468026E-0001 4.342437493468026E-0001
9.604912687080203E-0001
bi: 1.046562260264673E-0001 4.013974147759622E-0001 4.013974147759622E-0001
1.046562260264673E-0001
N = 4
xi: -8.611363115940526E-0001 -3.399810435848563E-0001 3.399810435848563E-0001
8.611363115940526E-0001
ai: 1.700536053357227E-0001 3.269491896014516E-0001 3.269491896014516E-0001
1.700536053357227E-0001
yi: -9.765602507375731E-0001 -6.402862174963100E-0001 0.000000000000000E+0000
6.402862174963100E-0001 9.765602507375731E-0001
bi: 6.297737366547301E-0002 2.667983404522844E-0001 3.464429818901364E-0001
2.667983404522844E-0001 6.297737366547301E-0002
N = 5
xi: -9.061798459386640E-0001 -5.384693101056831E-0001 0.000000000000000E+0000
5.384693101056831E-0001 9.061798459386640E-0001
ai: 1.152333166224734E-0001 2.410403392286476E-0001 2.829874178574912E-0001
2.410403392286476E-0001 1.152333166224734E-0001
yi: -9.840853600948425E-0001 -7.541667265708492E-0001 -2.796304131617832E-0001
2.796304131617832E-0001 7.541667265708492E-0001 9.840853600948425E-0001
bi: 4.258203675108183E-0002 1.868007965564927E-0001 2.728498019125589E-0001
2.728498019125589E-0001 1.868007965564927E-0001 4.258203675108183E-0002
N = 6
xi: -9.324695142031520E-0001 -6.612093864662645E-0001 -2.386191860831969E-0001
2.386191860831969E-0001 6.612093864662645E-0001 9.324695142031520E-0001
ai: 8.369444044690663E-0002 1.810719943231376E-0001 2.337708641169944E-0001
2.337708641169944E-0001 1.810719943231376E-0001 8.369444044690663E-0002
yi: -9.887032026126789E-0001 -8.213733408650279E-0001 -4.631182124753046E-0001
0.000000000000000E+0000 4.631182124753046E-0001 8.213733408650279E-0001
9.887032026126789E-0001
bi: 3.039615411981977E-0002 1.373206046344469E-0001 2.132096522719623E-0001
2.410725801734648E-0001 2.132096522719623E-0001 1.373206046344469E-0001
3.039615411981977E-0002
N = 7
xi: -9.491079123427585E-0001 -7.415311855993944E-0001 -4.058451513773972E-0001
0.000000000000000E+0000 4.058451513773972E-0001 7.415311855993944E-0001
9.491079123427585E-0001
ai: 6.309209262997855E-0002 1.406532597155259E-0001 1.903505780647854E-0001
2.094821410847278E-0001 1.903505780647854E-0001 1.406532597155259E-0001
6.309209262997855E-0002
yi: -9.914553711208126E-0001 -8.648644233597691E-0001 -5.860872354676911E-0001
-2.077849550078985E-0001 2.077849550078985E-0001 5.860872354676911E-0001
8.648644233597691E-0001 9.914553711208126E-0001
bi: 2.293532201052922E-0002 1.047900103222502E-0001 1.690047266392679E-0001
2.044329400752989E-0001 2.044329400752989E-0001 1.690047266392679E-0001
1.047900103222502E-0001 2.293532201052923E-0002
N = 8
xi: -9.602898564975362E-0001 -7.966664774136267E-0001 -5.255324099163290E-0001
-1.834346424956498E-0001 1.834346424956498E-0001 5.255324099163290E-0001
7.966664774136267E-0001 9.602898564975362E-0001
ai: 4.943939500213931E-0002 1.116463708268396E-0001 1.566526061681884E-0001
1.814000250680346E-0001 1.814000250680346E-0001 1.566526061681884E-0001
1.116463708268396E-0001 4.943939500213931E-0002
yi: -9.933798758817162E-0001 -8.941209068474564E-0001 -6.723540709451587E-0001
-3.607010979281320E-0001 0.000000000000000E+0000 3.607010979281320E-0001
6.723540709451587E-0001 8.941209068474564E-0001 9.933798758817162E-0001
bi: 1.782238332071036E-0002 8.248229893135833E-0002 1.362631092551722E-0001
1.720706085552113E-0001 1.844464057446916E-0001 1.720706085552113E-0001
1.362631092551722E-0001 8.248229893135833E-0002 1.782238332071036E-0002
N = 9
xi: -9.681602395076261E-0001 -8.360311073266358E-0001 -6.133714327005904E-0001
-3.242534234038089E-0001 0.000000000000000E+0000 3.242534234038089E-0001
6.133714327005904E-0001 8.360311073266358E-0001 9.681602395076261E-0001
ai: 3.963189516026125E-0002 9.079068168872639E-0002 1.300014068553412E-0001
1.564135277884839E-0001 1.648960128283494E-0001 1.564135277884839E-0001
1.300014068553412E-0001 9.079068168872639E-0002 3.963189516026126E-0002
yi: -9.946781606773402E-0001 -9.149635072496779E-0001 -7.344867651839338E-0001
-4.754624791124599E-0001 -1.642235636149868E-0001 1.642235636149868E-0001
4.754624791124599E-0001 7.344867651839338E-0001 9.149635072496779E-0001
9.946781606773402E-0001
bi: 1.430477564383894E-0002 6.651815594027414E-0002 1.117891346844183E-0001
1.452395883843662E-0001 1.628628274401151E-0001 1.628628274401151E-0001
1.452395883843662E-0001 1.117891346844183E-0001 6.651815594027414E-0002
1.430477564383894E-0002
N = 10
xi: -9.739065285171717E-0001 -8.650633666889845E-0001 -6.794095682990244E-0001
-4.333953941292472E-0001 -1.488743389816312E-0001 1.488743389816312E-0001
4.333953941292472E-0001 6.794095682990244E-0001 8.650633666889845E-0001
9.739065285171717E-0001
ai: 3.255816230796473E-0002 7.503967481091995E-0002 1.093871588022976E-0001
1.347092173114733E-0001 1.477391049013385E-0001 1.477391049013385E-0001
1.347092173114733E-0001 1.093871588022976E-0001 7.503967481091995E-0002
3.255816230796473E-0002
yi: -9.956571630258081E-0001 -9.301574913557082E-0001 -7.808177265864169E-0001
-5.627571346686047E-0001 -2.943928627014602E-0001 0.000000000000000E+0000
2.943928627014602E-0001 5.627571346686047E-0001 7.808177265864169E-0001
9.301574913557082E-0001 9.956571630258081E-0001
bi: 1.169463886737187E-0002 5.475589657435200E-0002 9.312545458369761E-0002
1.234919762620659E-0001 1.427759385770601E-0001 1.494455540029169E-0001
1.427759385770601E-0001 1.234919762620659E-0001 9.312545458369761E-0002
5.475589657435200E-0002 1.169463886737187E-0002
