(57)【特許請求の範囲】
【請求項１】
１次元以上の領域R上で関数F(・)により定義される力学系
xn+1 = F(xn)  (n≧1)
から生成された領域Rの要素の系列
x1，x2，…
がその極限密度分布関数
ρ(x)dx
についてエルゴード性を有するものに対して、
内積
<u,v> = ∫R u(x)v(x)ρ(x)dx
によりノルム演算
|v|2 = <v,v>
が定義されるヒルベルト空間Lの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルp  (p∈P)と各クラスの関数を指定するラベルq  (q∈Q(p))に対する基底関数
φp,q(・)
が、任意のx∈Rについて
φp,q(F(x)) = φp,q+1(x)
および
∫R φp,q(x)ρ(x)dx = 0
を満たすものを用いて、
関数A(・)の領域Rにおける積分
∫R A(x)dx
を求める数値積分装置であって、
系列x1，x2，…、ラベルp∈P，q∈Q(p)のそれぞれに対する係数Cp,q、関数A(・)の定義、ラベルp∈Pのそれぞれに対する係数Δp、関数B(・)の定義を記憶することができるＲＡＭ（Random Access Memory）、
前記ＲＡＭに、所定の初期値x1を記憶させ、前記ＲＡＭにすでに記憶済みのxj (j≧1)に関数F(・)を適用する漸化式
xj+1 = F(xj) (j≧1)
によりxj+1を計算して前記ＲＡＭに記憶させる生成部、
前記生成部による処理の後、前記ＲＡＭに記憶された系列x1，x2，…、関数A(・)の定義から、所定の整数Mを用いて、ラベルp∈P，q∈Q(p)のそれぞれについて、係数
Cp,q = (1/M)Σj=1M A(xj)φp,q(xj)/ρ(xi)
を計算して前記ＲＡＭに記憶させるC計算部、
前記C計算部による処理の後、前記ＲＡＭに記憶された係数Cp,qから、ラベルp∈Pのそれぞれについて、係数
Δp = -Σq∈Q(p) Cp,q
を計算して前記ＲＡＭに記憶させるΔ計算部、
前記Δ計算部による処理の後、終了条件が満たされない場合、前記ＲＡＭに記憶された係数Δp、関数A(・)の定義から、所定のq(P)∈Q(P)を用いて、関数B(・)の定義を
B(x) = A(x) + Σp∈P Δpφp,q(p)(x)ρ(x)
により計算して前記ＲＡＭに記憶させるB計算部、
前記B計算部による処理の後、前記ＲＡＭに記憶された関数B(・)の定義で前記ＲＡＭに記憶される関数A(・)の定義を置き換えて記憶させるA更新部、
前記C計算部による処理の後、前記ＲＡＭに記憶されたCp,qが前記終了条件を満たす場合、前記ＲＡＭに記憶された関数A(・)の定義から、所定の整数Nに対して
I = (1/N)Σi=1N A(xi)/ρ(xi)
を計算して、Iを数値積分結果として出力する積分部を備え、
前記終了条件は、前記ＲＡＭに記憶されたΔpについて、所定の誤差εに対して
Σp∈P Δp2 ≦ ε
が成立する場合に満たされ、
前記終了条件が満たされない場合、前記C計算部による処理を再度実行させる
ことを特徴とするもの。
【請求項２】
１次元以上の領域R上で関数F(・)により定義される力学系
xn+1 = F(xn) (n≧1)
から生成された領域Rの要素の系列
x1，x2，…
がその極限密度分布関数
ρ(x)dx
についてエルゴード性を有するものに対して、
内積
<u,v> = ∫R u(x)v(x)ρ(x)dx
によりノルム演算
|v|2 = <v,v>
が定義されるヒルベルト空間Lの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルp (p∈P)と各クラスの関数を指定するラベルq (q∈Q(p))に対する基底関数
φp,q(・)
が、任意のx∈Rについて
φp,q(F(x)) = φp,q+1(x)
および
∫R φp,q(x)ρ(x)dx = 0
を満たすものを用いて、
関数A(・)の領域Rにおける積分
∫R A(x)dx
を求める数値積分装置であって、
系列x1，x2，…、ラベルp∈P，q∈Q(p)のそれぞれに対する係数Cp,q、関数A(・)の定義、ラベルp∈Pのそれぞれに対する係数Δp、関数B(・)の定義を記憶することができるＲＡＭ（Random Access Memory）、
前記ＲＡＭに、所定の初期値x1を記憶させ、前記ＲＡＭにすでに記憶済みのxj (j≧1)に関数F(・)を適用する漸化式
xj+1 = F(xj) (j≧1)
によりxj+1を計算して前記ＲＡＭに記憶させる生成部、
前記生成部による処理の後、前記ＲＡＭに記憶された系列x1，x2，…、関数A(・)の定義から、所定の整数Mを用いて、ラベルp∈P，q∈Q(p)のそれぞれについて、係数
Cp,q = (1/M)Σj=1M A(xj)φp,q(xj)/ρ(xi)
を計算して前記ＲＡＭに記憶させるC計算部、
前記C計算部による処理の後、前記ＲＡＭに記憶された係数Cp,qから、ラベルp∈Pのそれぞれについて、係数
Δp = -Σq∈Q(p) Cp,q
を計算して前記ＲＡＭに記憶させるΔ計算部、
前記Δ計算部による処理の後、終了条件が満たされない場合、前記ＲＡＭに記憶された係数Δp、関数A(・)の定義から、所定のq(P)∈Q(P)を用いて、関数B(・)の定義を
B(x) = A(x) + Σp∈P Δpφp,q(p)(x)ρ(x)
により計算して前記ＲＡＭに記憶させるB計算部、
前記B計算部による処理の後、前記ＲＡＭに記憶された関数B(・)の定義で前記ＲＡＭに記憶される関数A(・)の定義を置き換えて記憶させるA更新部、
前記C計算部による処理の後、前記ＲＡＭに記憶されたCp,qが前記終了条件を満たす場合、前記ＲＡＭに記憶された関数A(・)の定義から、所定の整数Nに対して
I = (1/N)Σi=1N A(xi)/ρ(xi)
を計算して、Iを数値積分結果として出力する積分部を備え、
前記終了条件は、前記ＲＡＭに記憶されたΔpについて、所定の誤差εに対して
Σp∈P Δp2 ≦ ε
が成立する場合に満たされ、
前記終了条件が満たされない場合、前記C計算部による処理を再度実行させる
ことを特徴とするもの。
【請求項３】
請求項１または２に記載の数値積分装置であって、
前記領域RはD次元の領域であり、
所定のD個の正整数Pmax1，…，PmaxDについて
P =｛(p[1],…,p[D]) | 
p[1],…,p[D]は整数；
1≦p[1]≦Pmax1；…；
1≦p[D]≦PmaxD｝
であり、任意のp∈Pと、所定のD個の正整数Qmax1(p)，…，QmaxD(p)と、について
Q(p) =｛(q[1],…,q[D]) | 
q[1],…,q[D]は整数；
1≦q[1]≦Qmax1(p)；…；
1≦q[D]≦QmaxD(p)｝
である
ことを特徴とするもの。
【請求項４】
請求項１から３のいずれか１項に記載の数値積分装置であって、
前記力学系を定義する関数F(・)は、チェビシェフ多項式
T(a,cosθ) = cos(aθ) (a≧2)
により、領域Rの次元Dと、いずれかの整数a1，…，aD (a1≧2，…，aD≧2)と、y=(y1，…，yD)∈Rとについて、
F(y) = (T(a1,y1)，…，T(aD,yD))
と定義され、
前記極限密度分布関数ρ(・)は、
ρ(x) = 1/〔πD(1-y12)1/2…(1-yD2)1/2〕
と定義され、
p = (p[1]，…，p[D])；
q = (q[1]，…，q[D])
に対して、前記基底関数φp,q(・)は、
φp,q(y) = 2D/2T(p[1]a1q[1],y1)…T(p[D]aDq[D],yD)
と定義される
ことを特徴とするもの。
【請求項５】
請求項１から４のいずれか１項に記載の数値積分装置であって、
前記C計算部による処理を再度実行させるのに先立って、前記生成部に、前記所定の初期値とは異なる値を初期値として与えて系列を新たに生成させて前記ＲＡＭに記憶させ、前記C計算部による処理を再度実行させる際の計算では当該新たに生成されて前記ＲＡＭに記憶された系列を用いる
ことを特徴とするもの。
【請求項６】
請求項１から５のいずれか１項に記載の数値積分装置であって、
前記積分部における計算に先立って、前記生成部に、前記所定の初期値とは異なる値を初期値として与えて系列を新たに生成させて前記ＲＡＭに記憶させ、前記積分部における計算では当該新たに生成されて前記ＲＡＭに記憶された系列を用いる
ことを特徴とするもの。
【請求項７】
１次元以上の領域R上で関数F(・)により定義される力学系
xn+1 = F(xn) (n≧1)
から生成された領域Rの要素の系列
x1，x2，…
がその極限密度分布関数
ρ(x)dx
についてエルゴード性を有するものに対して、
内積
<u,v> = ∫R u(x)v(x)ρ(x)dx
によりノルム演算
|v|2 = <v,v>
が定義されるヒルベルト空間Lの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルp (p∈P)と各クラスの関数を指定するラベルq (q∈Q(p))に対する基底関数
φp,q(・)
が、任意のx∈Rについて
φp,q(F(x)) = φp,q+1(x)
および
∫R φp,q(x)ρ(x)dx = 0
を満たすものを用いて、
関数A(・)の領域Rにおける積分
∫R A(x)dx
を求める数値積分方法であって、系列x1，x2，…、ラベルp∈P，q∈Q(p)のそれぞれに対する係数Cp,q、関数A(・)の定義、ラベルp∈Pのそれぞれに対する係数Δp、関数B(・)の定義を記憶することができるＲＡＭ（Random Access Memory）、生成部、C計算部、Δ計算部、B計算部、A更新部、積分部を有する数値積分装置にて実行され、
前記生成部が、前記ＲＡＭに、所定の初期値x1を記憶させ、前記ＲＡＭにすでに記憶済みのxj (j≧1)に関数F(・)を適用する漸化式
xj+1 = F(xj) (j≧1)
によりxj+1を計算して前記ＲＡＭに記憶させる生成工程、
前記生成工程の後、前記C計算部が、前記ＲＡＭに記憶された系列x1，x2，…、関数A(・)の定義から、所定の整数Mを用いて、ラベルp∈P，q∈Q(p)のそれぞれについて、係数
Cp,q = (1/M)Σj=1M A(xj)φp,q(xj)/ρ(xi)
を計算して前記ＲＡＭに記憶させるC計算工程、
前記C計算工程の後、前記Δ計算部が、前記ＲＡＭに記憶された係数Cp,qから、ラベルp∈Pのそれぞれについて、係数
Δp = -Σq∈Q(p) Cp,q
を計算して前記ＲＡＭに記憶させるΔ計算工程、
前記Δ計算工程の後、前記B計算部が、終了条件が満たされない場合、前記ＲＡＭに記憶された係数Δp、関数A(・)の定義から、所定のq(P)∈Q(P)を用いて、関数B(・)の定義を
B(x) = A(x) + Σp∈P Δpφp,q(p)(x)ρ(x)
により計算して前記ＲＡＭに記憶させるB計算工程、
前記B計算工程の後、前記A更新部が、前記ＲＡＭに記憶された関数B(・)の定義で前記ＲＡＭに記憶される関数A(・)の定義を置き換えて記憶させるA更新工程、
前記C計算工程の後、前記積分部が、前記ＲＡＭに記憶されたCp,qが当該終了条件を満たす場合、前記ＲＡＭに記憶された関数A(・)の定義から、所定の整数Nに対して
I = (1/N)Σi=1N A(xi)/ρ(xi)
を計算して、Iを数値積分結果として出力する積分工程を備え、
前記終了条件は、前記ＲＡＭに記憶されたΔpについて、所定の誤差εに対して
Σp∈P Δp2 ≦ ε
が成立する場合に満たされ、
前記終了条件が満たされない場合、前記C計算工程を再度実行する
ことを特徴とする方法。
【請求項８】
１次元以上の領域R上で関数F(・)により定義される力学系
xn+1 = F(xn) (n≧1)
から生成された領域Rの要素の系列
x1，x2，…
がその極限密度分布関数
ρ(x)dx
についてエルゴード性を有するものに対して、
内積
<u,v> = ∫R u(x)v(x)ρ(x)dx
によりノルム演算
|v|2 = <v,v>
が定義されるヒルベルト空間Lの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルp (p∈P)と各クラスの関数を指定するラベルq (q∈Q(p))に対する基底関数
φp,q(・)
が、任意のx∈Rについて
φp,q(F(x)) = φp,q+1(x)
および
∫R φp,q(x)ρ(x)dx = 0
を満たすものを用いて、
関数A(・)の領域Rにおける積分
∫R A(x)dx
を求める数値積分方法であって、系列x1，x2，…、ラベルp∈P，q∈Q(p)のそれぞれに対する係数Cp,q、関数A(・)の定義、ラベルp∈Pのそれぞれに対する係数Δp、関数B(・)の定義を記憶することができるＲＡＭ（Random Access Memory）、生成部、C計算部、Δ計算部、B計算部、A更新部、積分部を有する数値積分装置にて実行され、
前記生成部が、前記ＲＡＭに、所定の初期値x1を記憶させ、前記ＲＡＭにすでに記憶済みのxj (j≧1)に関数F(・)を適用する漸化式
xj+1 = F(xj) (j≧1)
によりxj+1を計算して前記ＲＡＭに記憶させる生成工程、
前記生成工程の後、前記C計算部が、前記ＲＡＭに記憶された系列x1，x2，…、関数A(・)の定義から、所定の整数Mを用いて、ラベルp∈P，q∈Q(p)のそれぞれについて、係数
Cp,q = (1/M)Σj=1M A(xj)φp,q(xj)/ρ(xi)
を計算して前記ＲＡＭに記憶させるC計算工程、
前記C計算工程の後、前記Δ計算部が、前記ＲＡＭに記憶された係数Cp,qから、ラベルp∈Pのそれぞれについて、係数
Δp = -Σq∈Q(p) Cp,q
を計算して前記ＲＡＭに記憶させるΔ計算工程、
前記Δ計算工程の後、前記B計算部が、終了条件が満たされない場合、前記ＲＡＭに記憶された係数Δp、関数A(・)の定義から、所定のq(P)∈Q(P)を用いて、関数B(・)の定義を
B(x) = A(x) + Σp∈P Δpφp,q(p)(x)ρ(x)
により計算して前記ＲＡＭに記憶させるB計算工程、
前記B計算工程の後、前記A更新部が、前記ＲＡＭに記憶された関数B(・)の定義で前記ＲＡＭに記憶される関数A(・)の定義を置き換えて記憶させるA更新工程、
前記C計算工程の後、前記積分部が、前記ＲＡＭに記憶されたCp,qが当該終了条件を満たす場合、前記ＲＡＭに記憶された関数A(・)の定義から、所定の整数Nに対して
I = (1/N)Σi=1N A(xi)/ρ(xi)
を計算して、Iを数値積分結果として出力する積分工程を備え、
前記終了条件は、前記ＲＡＭに記憶されたΔpについて、所定の誤差εに対して
Σp∈P Δp2 ≦ ε
が成立する場合に満たされ、
前記終了条件が満たされない場合、前記C計算工程を再度実行する
ことを特徴とする方法。
【請求項９】
請求項７または８に記載の数値積分方法であって、
前記領域RはD次元の領域であり、
所定のD個の正整数Pmax1，…，PmaxDについて
P =｛(p[1],…,p[D]) | 
p[1],…,p[D]は整数；
1≦p[1]≦Pmax1；…；
1≦p[D]≦PmaxD｝
であり、任意のp∈Pと、所定のD個の正整数Qmax1(p)，…，QmaxD(p)と、について
Q(p) =｛(q[1],…,q[D]) | 
q[1],…,q[D]は整数；
1≦q[1]≦Qmax1(p)；…；
1≦q[D]≦QmaxD(p)｝
である
ことを特徴とする方法。
【請求項１０】
請求項７から９のいずれか１項に記載の数値積分方法であって、
前記力学系を定義する関数F(・)は、チェビシェフ多項式
T(a,cosθ) = cos(aθ) (a≧2)
により、領域Rの次元Dと、いずれかの整数a1，…，aD (a1≧2，…，aD≧2)と、y=(y1，…，yD)∈Rとについて、
F(y) = (T(a1,y1)，…，T(aD,yD))
と定義され、
前記極限密度分布関数ρ(・)は、
ρ(x) = 1/〔πD(1-y12)1/2…(1-yD2)1/2〕
と定義され、
p = (p[1]，…，p[D])；
q = (q[1]，…，q[D])
に対して、前記基底関数φp,q(・)は、
φp,q(y) = 2D/2T(p[1]a1q[1],y1)…T(p[D]aDq[D],yD)
と定義される
ことを特徴とする方法。
【請求項１１】
請求項７から１０のいずれか１項に記載の数値積分方法であって、
前記C計算工程を再度実行させるのに先立って、前記生成部に、前記所定の初期値とは異なる値を初期値として与えて系列を新たに生成させて前記ＲＡＭに記憶させ、前記C計算工程を再度実行させる際の計算では当該新たに生成されて前記ＲＡＭに記憶された系列を用いる
ことを特徴とするもの。
【請求項１２】
請求項７から１１のいずれか１項に記載の数値積分方法であって、
前記積分工程における計算に先立って、前記生成部に、前記所定の初期値x1とは異なる値を初期値として与えて系列を新たに生成させて前記ＲＡＭに記憶させ、前記積分工程における計算では当該新たに生成されて前記ＲＡＭに記憶された系列を用いる
ことを特徴とする方法。
【請求項１３】
コンピュータを請求項１から６のいずれか１項に記載の数値積分装置として機能させることを特徴とするプログラム。
【請求項１４】
コンピュータに、請求項７から１２のいずれか１項に記載の数値積分方法を実行させることを特徴とするプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、数値積分装置、数値積分方法、ならびに、これらをコンピュータ上で実現するプログラムに関する。
【０００２】
【従来の技術】
従来から、乱数を用いて数値積分を行う方法としてモンテカルロ法という手法が知られている。特に多次元の数値積分の分野においては、事実上モンテカルロ法に頼らなければ数値積分ができないことが多い。
【０００３】
モンテカルロ法では、乱数の利用を前提とするため、性質が良好な乱数が必要とされる。発明者らは、このような乱数としてカオス乱数を利用する技術について研究を進めている。
【０００４】
カオス乱数では、最も基本的な乱数である線形合同法と同様に、漸化式を用いて乱数列を発生させる。すなわち、
ｘｉ＋１＝Ｆ（ｘｉ）  （ｉ≧０）
なる漸化式を用いて得た数列
ｘ０，ｘ１，ｘ２，…
を擬似乱数列として扱うのである。この漸化式は、一般に「力学系」と呼ばれる。
【０００５】
カオス乱数では、漸化式に用いられる写像（関数）Ｆとしては、
ｆ（ｘ）＝４ｘ（１−ｘ）
のように定義されるウラム＝フォン・ノイマン写像、
ｆ（ｘ）＝３（３−４ｘ）２
のように定義されるキュービック写像、
ｆ（ｘ）＝ｘ（５−２０ｘ＋１６ｘ２）２
のように定義されるクインティック写像などが知られている。
【０００６】
これらは、
Ｕ（ａ，ｓｉｎ２θ）＝ｓｉｎ２（ａθ）
のように定義される多項式Ｕ（ａ，・）の一種である。
【０００７】
このほか、
Ｔ（ａ，ｃｏｓθ）＝ｃｏｓ（ａθ）
のように定義されるチェビシェフ多項式Ｔ（ａ，・）も利用できる。
【０００８】
このような写像を利用して生成されるカオス乱数には、以下のような性質があることが判明しており、その重要性はますます高まってきている。
【０００９】
（１）出力される乱数列には周期がないため、繰り返し同じ列が出力されることがない。
（２）乱数の種（漸化式に与える初期値）として有理数を与えると、得られる乱数列に含まれる数がいずれも有理数になる。
（３）乱数の分布を表す密度関数が既知の解析的関数である。
【００１０】
生成される乱数の分布は、具体的には、
ρ（ｘ）＝１／〔π（ｘ（１−ｘ））１／２〕
のように表現される。
【００１１】
また、パラメータを有する有理写像も漸化式として用いることができ、このような有理写像としてカツラ＝フクダ写像、一般化ウラム＝フォン・ノイマン写像、一般化キュービック写像、一般化チェビシェフ写像などがある。このうち、一般化ウラム＝フォン・ノイマン写像は、以下のように表現することができる。
ｆ（ｎ，ｍ，ｘ）＝４ｘ（１−ｘ）（１−ｎｘ）（１−ｍｘ）／（１−Ａｘ＋Ｂｘ２＋Ｃｘ４）
ただし
Ａ＝２（ｎ＋ｍ＋ｎｍ）
Ｂ＝８ｎｍ
Ｃ＝ｎ２＋ｍ２−２ｎｍ−２ｎ２ｍ−２ｎｍ２＋ｎ２ｍ２
【００１２】
この一般化ウラム＝フォン・ノイマン写像を用いて上記の漸化式により乱数列を得た場合、その分布は、同じパラメータを有する
ρ（ｘ）＝１／〔Ｋ（ｎ，ｍ）（ｘ（１−ｘ）（１−ｋｘ）（１−ｍｘ））１／２〕
ただし
Ｋ（ｎ，ｍ）＝∫０１ｄｕ／〔（１−ｕ２）（１−ｎｕ２）（１−ｍｕ２）〕１／２
なる密度関数で表現される。
【００１３】
なお、カツラ＝フクダ写像は、一般化ウラム＝フォン・ノイマン写像において、ｍ＝０とおいたものである。
【００１４】
これらの有理写像を用いた乱数を生成する手法については、本願の発明者らによる出願に係る特開平１０−２８３３４４号公報に開示されている。
【００１５】
【発明が解決しようとする課題】
さて、このような乱数を用いるモンテカルロ法による積分計算では、積分値の誤差εの分散Ｖ（ε）は、発生させた乱数の個数Ｎに反比例するのが一般的である。したがって、一般に積分計算の収束は低速であると言わざるを得ない。
【００１６】
そこで、モンテカルロ法による積分計算の収束が高速となるような技術が強く望まれている。
【００１７】
本発明は、以上のような問題を解決するためになされたもので、数値積分装置、数値積分方法、ならびに、コンピュータ上でこれらを実現するプログラムを提供することを目的とする。
【００１８】
【課題を解決するための手段】
以上の目的を達成するため、本発明の原理にしたがって、下記の発明を開示する。
【００１９】
本発明の第１の観点にかかる数値積分装置は、１次元以上の領域Ｒ上で関数Ｆ（・）により定義される力学系
ｘｎ＋１＝Ｆ（ｘｎ）  （ｎ≧１）
から生成された領域Ｒの要素の系列
ｘ１，ｘ２，…
がその極限密度分布関数
ρ（ｘ）ｄｘ
についてエルゴード性を有するものに対して、
【００２０】
内積
＜ｕ，ｖ＞＝∫Ｒｕ（ｘ）ｖ（ｘ）ρ（ｘ）ｄｘ
によりノルム演算
｜ｖ｜２＝＜ｖ，ｖ＞
が定義されるヒルベルト空間Ｌの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルｐ（ｐ∈Ｐ）と各クラスの関数を指定するラベルｑ（ｑ∈Ｑ（ｐ））に対する基底関数
φｐ，ｑ（・）
が、任意のｘ∈Ｒについて
φｐ，ｑ（Ｆ（ｘ））＝φｐ，ｑ＋１（ｘ）
および
∫Ｒφｐ，ｑ（ｘ）ρ（ｘ）ｄｘ＝０
を満たすものを用いて、
【００２１】
関数Ａ（・）の領域Ｒにおける積分
∫ＲＡ（ｘ）ｄｘ
を求める。
【００２２】
そして、生成部と、Ｃ計算部と、Δ計算部と、Ｂ計算部と、Ａ更新部と、積分部と、を備え、以下のように構成する。
【００２３】
すなわち、生成部は、所定の初期値ｘ１から、漸化式
ｘｊ＋１＝Ｆ（ｘｊ）  （ｊ≧１）
により系列
ｘ１，ｘ２，…
を生成する。
【００２４】
一方、Ｃ計算部は、所定の整数Ｍを用いて、ラベルｐ∈Ｐ，ｑ∈Ｑ（ｐ）のそれぞれについて、
Ｃｐ，ｑ＝（１／Ｍ）Σｊ＝１ＭＡ（ｘｊ）φｐ，ｑ（ｘｊ）／ρ（ｘｉ）
を計算する。
【００２５】
さらに、Δ計算部は、計算されたＣｐ，ｑが終了条件を満たさない場合、
Δｐ＝−Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ
を計算する。
【００２６】
そして、Ｂ計算部は、所定のｑ（Ｐ）∈Ｑ（Ｐ）を用いて、関数
Ｂ（ｘ）＝Ａ（ｘ）＋Σｐ∈ＰΔｐφｐ，ｑ（ｐ）（ｘ）ρ（ｘ）
を計算する。
【００２７】
さらに、Ａ更新部は、計算された関数Ｂ（ｘ）で関数Ａ（ｘ）を置き換える。
【００２８】
一方、積分部は、計算されたＣｐ，ｑが終了条件を満たす場合、所定の整数Ｎに対して
Ｉ＝（１／Ｎ）Σｉ＝ｌＮＡ（ｘｉ）／ρ（ｘｉ）
を計算し、Ｉを数値積分結果とする。
【００２９】
また、本発明の数値積分装置は、以下のように構成することができる。
【００３０】
すなわち、Δ計算部は、終了条件を満たすか否かにかかわらずΔｐを計算する。
【００３１】
一方、終了条件は、計算されたΔｐについて、所定の誤差εに対して
Σｐ∈ＰΔｐ２≦ε
が成立する場合に満たされる。
【００３２】
さらに、Ｂ計算部は、終了条件が満たされない場合、Ｂ（ｘ）を計算する。
【００３３】
また、本発明の数値積分装置は、以下のように構成することができる。
【００３４】
すなわち、Δ計算部は、終了条件を満たすか否かにかかわらずΔｐを計算する。
【００３５】
一方、終了条件は、任意のラベルｐ∈Ｐについて、所定の誤差εに対して、
｜Δｐ｜≦ε
が成立する場合に満たされる。
【００３６】
さらに、Ｂ計算部は、終了条件が満たされない場合、Ｂ（ｘ）を計算する。
【００３７】
また、本発明の数値積分装置は、以下のように構成することができる。
【００３８】
すなわち、領域ＲはＤ次元の領域である。
【００３９】
一方、所定のＤ個の正整数Ｐｍａｘｌ，…，ＰｍａｘＤについて
Ｐ＝｛（ｐ［１］，…，ｐ［Ｄ］）｜
ｐ［１］，…，ｐ［Ｄ］は整数；
１≦ｐ［１］≦Ｐｍａｘｌ…；
１≦ｐ［Ｄ］≦ＰｍａｘＤ｝
であり、任意のｐ∈Ｐと、所定のＤ個の正整数Ｑｍａｘ１（ｐ），…，ＱｍａｘＤ（ｐ）と、について
Ｑ（ｐ）＝｛（ｑ［１］，…，ｑ［Ｄ］）｜
ｑ［１］，…，ｑ［Ｄ］は整数；
１≦ｑ［１］≦Ｑｍａｘｌ（ｐ）；…；
１≦ｑ［Ｄ］≦ＱｍａｘＤ（ｐ）｝
である。
【００４０】
また、本発明の数値積分装置は、以下のように構成することができる。
【００４１】
すなわち、力学系を定義する関数Ｆ（・）は、チェビシェフ多項式
Ｔ（ａ，ｃｏｓθ）＝ｃｏｓ（ａθ）  （ａ≧２）
により、領域Ｒの次元Ｄと、いずれかの整数ａ１，…，ａＤ（ａ１≧２，…，ａＤ≧２）と、ｙ＝（ｙ１，…，ｙＤ）∈Ｒとについて、
Ｆ（ｙ）＝（Ｔ（ａ１，ｙ１），…，Ｔ（ａＤ，ｙＤ））
と定義される。
【００４２】
一方、極限密度分布関数ρ（・）は、
ρ（ｘ）＝１／〔πＤ（１−ｙｌ２）１／２…（１−ｙＤ２）１／２〕
と定義される。
【００４３】
さらに、
ｐ＝（ｐ［１］，…，ｐ［Ｄ］）；
ｑ＝（ｑ［１］…ｑ［Ｄ］）
に対して、基底関数φｐ，ｑ（・）は、
φｐ，ｑ（ｙ）＝２Ｄ／２Ｔ（ｐ［１］ａｌｑ［１］，ｙｌ）…Ｔ（ｐ［Ｄ］ａＤｑ［Ｄ］，ｙＤ）
と定義される。
【００４４】
また、本発明の数値積分装置は、初期値付与部と、平均出力部と、をさらに備え、以下のように構成することができる。
【００４５】
すなわち、初期値付与部は、生成部に、互いに異なる初期値を当該所定の初期値ｘｌとして与える。
【００４６】
一方、平均出力部は、互いに異なる初期値のそれぞれを与えた結果積分部により計算されるＩの相加平均を、数値積分結果として出力する。
【００４７】
また、本発明の数値積分装置において終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算部により計算されるΔｐの自乗和の総和が、所定の誤差ε以下である場合に満たされるように構成することができる。
【００４８】
また、本発明の数値積分装置において終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算部により計算されるΔｐのいずれもが、所定の誤差ε以下である場合に満たされるように構成することができる。
【００４９】
また、本発明の数値積分装置において互いに異なる初期値のそれぞれに対して、生成部、Ｃ計算部、Δ計算部、Ｂ計算部、Ａ更新部、ならびに、積分部により行われる処理は、並列に行われるように構成することができる。
【００５０】
本発明の他の観点に係る数値積分方法は、１次元以上の領域Ｒ上で関数Ｆ（・）により定義される力学系
ｘｎ＋１＝Ｆ（ｘｎ）  （ｎ≧１）
から生成された領域Ｒの要素の系列
ｘ１，ｘ２…
がその極限密度分布関数
ρ（ｘ）ｄｘ
についてエルゴード性を有するものに対して、
【００５１】
内積
＜ｕ，ｖ＞＝∫Ｒｕ（ｘ）ｖ（ｘ）ρ（ｘ）ｄｘ
によりノルム演算
｜ｖ｜２＝＜ｖ，ｖ＞
が定義されるヒルベルト空間Ｌの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルｐ（ｐ∈Ｐ）と各クラスの関数を指定するラベルｑ（ｑ∈Ｑ（ｐ））に対する基底関数
φｐ，ｑ（・）
が、任意のｘ∈Ｒについて
φｐ，ｑ（Ｆ（ｘ））＝φｐ，ｑ＋１（ｘ）
および
∫Ｒφｐ，ｑ（ｘ）ρ（ｘ）ｄｘ＝０
を満たすものを用いて、
【００５２】
関数Ａ（・）の領域Ｒにおける積分
∫ＲＡ（ｘ）ｄｘ
を求める。
【００５３】
そして、生成工程と、Ｃ計算工程と、Δ計算工程と、Ｂ計算工程と、Ａ更新工程と、積分工程と、を備え、以下のように構成する。
【００５４】
すなわち、生成工程では、所定の初期値ｘ１から、漸化式
ｘｊ＋１＝Ｆ（ｘｊ）  （ｊ≧１）
により系列
ｘ１，ｘ２，…
を生成する。
【００５５】
一方、Ｃ計算工程では、所定の整数Ｍを用いて、ラベルｐ∈Ｐ，ｑ∈Ｑ（ｐ）のそれぞれについて、
Ｃｐ，ｑ＝（１／Ｍ）Σｊ＝ｌＭＡ（ｘｊ）φｐ，ｑ（ｘｊ）／ρ（ｘｉ）
を計算する。
【００５６】
さらに、Δ計算工程では、計算されたＣｐ，ｑが終了条件を満たさない場合、
Δｐ＝−Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ
を計算する。
【００５７】
そして、Ｂ計算工程では、所定のｑ（Ｐ）∈Ｑ（Ｐ）を用いて、関数
Ｂ（ｘ）＝Ａ（ｘ）＋Σｐ∈ＰΔｐφｐ，ｑ（Ｐ）（ｘ）ρ（ｘ）
を計算する。
【００５８】
さらに、Ａ更新工程では、計算された関数Ｂ（ｘ）で関数Ａ（ｘ）を置き換える。
【００５９】
一方、積分工程では、計算されたＣｐ，ｑが終了条件を満たす場合、所定の整数Ｎに対して
Ｉ＝（１／Ｎ）Σｉ＝ｌＮＡ（ｘｉ）／ρ（ｘｉ）
を計算し、Ｉを数値積分結果とする。
【００６０】
また、本発明の数値積分方法は、以下のように構成することができる。
【００６１】
すなわち、Δ計算工程では、終了条件を満たすか否かにかかわらずΔｐを計算する。
【００６２】
一方、終了条件は、計算されたΔｐについて、所定の誤差εに対して
Σｐ∈ＰΔｐ２≦ε
が成立する場合に満たされる。
【００６３】
さらに、Ｂ計算工程では、終了条件が満たされない場合、Ｂ（ｘ）を計算する。
【００６４】
また、本発明の数値積分方法は、以下のように構成することができる。
【００６５】
すなわち、Δ計算工程では、終了条件を満たすか否かにかかわらずΔｐを計算する。
【００６６】
一方、終了条件は、任意のラベルｐ∈Ｐについて、所定の誤差εに対して、
｜Δｐ｜≦ε
が成立する場合に満たされる。
【００６７】
さらに、Ｂ計算工程では、終了条件が満たされない場合、Ｂ（ｘ）を計算する。
【００６８】
また、本発明の数値積分方法は、以下のように構成することができる。
【００６９】
すなわち、領域ＲはＤ次元の領域である。
【００７０】
一方、所定のＤ個の正整数Ｐｍａｘｌ，…，ＰｍａｘＤについて
Ｐ＝｛（ｐ［１］，…，ｐ［Ｄ］）｜
ｐ［１］，…，ｐ［Ｄ］は整数；
１≦ｐ［１］≦Ｐｍａｘｌ；…；
１≦ｐ［Ｄ］≦ＰｍａｘＤ｝
であり、任意のｐ∈Ｐと、所定のＤ個の正整数Ｑｍａｘｌ（ｐ），…，ＱｍａｘＤ（ｐ）と、について
Ｑ（ｐ）＝｛（ｑ［１］，…，ｑ［Ｄ］）｜
ｑ［１］，…，ｑ［Ｄ］は整数；
１≦ｑ［１」≦Ｑｍａｘｌ（ｐ）；…；
１≦ｑ［Ｄ」≦ＱｍａｘＤ（ｐ）｝
である。
【００７１】
また、本発明の数値積分方法は、以下のように構成することができる。
【００７２】
すなわち、力学系を定義する関数Ｆ（・）は、チェビシェフ多項式
Ｔ（ａ，ｃｏｓθ）＝ｃｏｓ（ａθ）  （ａ≧２）
により、領域Ｒの次元Ｄと、いずれかの整数ａｌ，…，ａＤ（ａｌ≧２，…，ａＤ≧２）と、ｙ＝（ｙｌ，…，ｙＤ）∈Ｒとについて、
Ｆ（ｙ）＝（Ｔ（ａｌ，ｙｌ），…，Ｔ（ａＤ，ｙＤ））
と定義される。
【００７３】
一方、極限密度分布関数ρ（・）は、
ρ（ｘ）＝１／〔πＤ（１−ｙｌ２）１／２…（１−ｙＤ２）１／２〕
と定義される。
【００７４】
さらに、
ｐ＝（ｐ［１］，…，ｐ［Ｄ］）；
ｑ＝（ｑ［１］，…，ｑ［Ｄ］）
に対して、基底関数φｐ，ｑ（・）は、
φｐ，ｑ（ｙ）＝２Ｄ／２Ｔ（ｐ［１］ａｌｑ［１］，ｙｌ）…Ｔ（ｐ［Ｄ］ａＤｑ［Ｄ］，ｙＤ）
と定義される。
【００７５】
また、本発明の数値積分方法は、初期値付与工程と、平均出力工程と、をさらに備え、以下のように構成することができる。
【００７６】
すなわち、初期値付与工程では、生成工程に、互いに異なる初期値を当該所定の初期値ｘｌとして与える。
【００７７】
一方、平均出力工程では、互いに異なる初期値のそれぞれを与えた結果積分工程にて計算されるＩの相加平均を、数値積分結果として出力する。
【００７８】
また、本発明の数値積分方法において、終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算工程にて計算されるΔｐの自乗和の総和が、所定の誤差ε以下である場合に満たされるように構成することができる。
【００７９】
また、本発明の数値積分方法において、終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算工程にて計算されるΔｐのいずれもが、所定の誤差ε以下である場合に満たされるように構成することができる。
【００８０】
また、本発明の数値積分方法において、互いに異なる初期値のそれぞれに対して、生成工程、Ｃ計算工程、Δ計算工程、Ｂ計算工程、Ａ更新工程、ならびに、積分工程にて行われる処理は、並列に行われるように構成することができる。
【００８１】
本発明の他の観点に係るプログラムは、コンピュータ（ＤＳＰ  （Ｄｉｇｉｔａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｏｒ）、ＦＰＧＡ（Ｆｉｅｌｄ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）を含む。）を上記の数値積分装置として機能させ、もしくは、コンピュータに上記の数値積分方法を実行させるように構成する。
【００８２】
当該プログラムは、コンパクトディスク、フレキシブルディスク、ハードディスク、光磁気ディスク、ディジタルビデオディスク、磁気テープ、半導体メモリなどのコンピュータ読取可能な情報記録媒体に記録することができる。
【００８３】
また、当該プログラムは、コンピュータ通信網を介して、これが実行されるコンピュータとは独立に配布・販売することができる。
【００８４】
また、当該プログラムを記録したコンピュータ読取可能な情報記録媒体は、これが実行されるコンピュータとは独立に配布・販売することができる。
【００８５】
【発明の実施の形態】
以下に本発明の一実施形態を説明する。なお、以下に説明する実施形態は説明のためのものであり、本願発明の範囲を制限するものではない。したがって、当業者であればこれらの各要素もしくは全要素をこれと均等なものに置換した実施形態を採用することが可能であるが、これらの実施形態も本願発明の範囲に含まれる。
【００８６】
（理論的背景）
上記のように、漸化式で示される力学系によって発生される乱数を用いてモンテカルロ法による関数Ａ（・）の領域Ｒにおける数値積分を行う場合、当該乱数
ｘ０，ｘ１，ｘ２，…
の極限分布がρ（ｘ）のように表現できるとすると、かりに
∫ＲＳ（ｘ）ρ（ｘ）ｄｘ＝０
を満たす関数Ｓ（ｘ）が見つかれば、
∫ＲＡ（ｘ）ｄｘ
＝∫Ｒ〔Ａ（ｘ）／ρ（ｘ）＋Ｓ（ｘ）〕ρ（ｘ）ｄｘ
≒（１／Ｎ）Σｉ＝０Ｎ−１〔Ａ（ｘｉ）／ρ（ｘｉ）＋Ｓ（ｘｉ）〕ρ（ｘｉ）〕
が成立する。最後の≒は、モンテカルロ法による近似を意味する。
【００８７】
さて、
∫Ｒ｜Ａ（ｘ）｜２ｄｘ＜∞
であれば、すなわち、左辺の積分が有界であれば、上記のモンテカルロ法によって得られた数値積分値の誤差の分散Ｖ（ε）が、
Ｖ（ε）∝１／Ｎ２
となるような関数Ｓ（・）が必ず存在することが、発明者らによって証明されている。
【００８８】
したがって、繰り返し計算において、上記の条件を満たしつつ、収束が早くなるように、関数Ｓ（・）を構成すればよいこととなる。
【００８９】
ここで、内積
＜ｕ，ｖ＞＝∫Ｒｕ（ｘ）ｖ（ｘ）ρ（ｘ）ｄｘ
によりノルム演算
｜ｖ｜２＝＜ｖ，ｖ＞
が定義されるヒルベルト空間Ｌの正規直交基底（ルベークスペクトラム）であって、当該正規直交基底のクラスを指定するラベルｐ（ｐ∈Ｐ）と各クラスの関数を指定するラベルｑ（ｑ∈Ｑ（ｐ））に対する基底関数
φｐ，ｑ（・）
が、任意のｘ∈Ｒについて
φｐ，ｑ（Ｆ（ｘ））＝φｐ，ｑ＋１（ｘ）
および
∫Ｒφｐ，ｑ（ｘ）ρ（ｘ）ｄｘ＝０
を満たすものを用いることを考える。
【００９０】
まず、これは正規直交基底であるから、
∫Ｒφｐ，ｉ（ｘ）φｐ，ｊ（ｘ）ρ（ｘ）ｄｘ＝δｉ，ｊ
が成立する。ただし、δｉ，ｊはクロネッカーのデルタ関数である。
【００９１】
また、関数Ａ（・）は、
Ａ（ｘ）＝Σｐ∈ＰΣｑ∈Ｑ（ｐ）Ｃｐ，ｑφｐ，ｑ（ｘ）
のように展開できる。
【００９２】
尚、実際には、展開を解析的に行うのではなく、こちらでも乱数を用いて、
Ｃｐ，ｑ＝（１／Ｍ）Σｊ＝ｌＭＡ（ｘｊ）φｐ，ｑ（ｘｊ）／ρ（ｘｉ）
のように展開することとなる。
【００９３】
これらと上記の性質とを比較すると、すべてのｐについて
Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ＝０
になれば、
Ｖ（ε）∝１／Ｎ２
が成立する。
【００９４】
このφｐ，ｑ（・）やρ（ｘ）に相当する関数はすでに存在する。たとえば、後者については、チェビシェフ多項式により生成されるカオス乱数を用いればよい。
【００９５】
上記のように、非積分関数をＡ（ｘ）からＡ（ｘ）＋Ｓ（ｘ）ρ（ｘ）に修正することによって、収束が高速になるのであるから、Ａ（ｘ）を展開した結果得られた係数Ｃｐ，ｑを修正することによって、収束を高速にすることができる。
【００９６】
実際にＡ（ｘ）を展開してみると、
Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ＝０
は成立しない。そこで、
Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ＋Δｐ＝０
となるようなΔｐを考え、所定のｑ（Ｐ）∈Ｑ（Ｐ）を用いて、以下のような関数Ｂ（・）について検討することとする。
Ｂ（ｘ）＝Ａ（ｘ）＋Σｐ∈ＰΔｐφｐ，ｑ（ｐ）（ｘ）ρ（ｘ）
【００９７】
さて、上述したような直交性の観点からすると、Ａ（ｘ）をモンテカルロ法で数値積分することと、Ｂ（ｘ）をモンテカルロ法で数値積分することとは、等価であると考えられる。すなわち、同じ結果が得られる。
【００９８】
しかも、Ｂ（ｘ）を
Ｂ（ｘ）＝Σｐ∈ＰΣｑ∈Ｑ（ｐ）Ｃ’ｐ，ｑφｐ，ｑ（ｘ）
のように展開した場合、各係数について
Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ
Σｑ∈Ｑ（ｐ）Ｃ’ｐ，ｑ
の２つを比較すると、後者の方がより０に近いことが期待される。したがって、Ａ（ｘ）をＢ（ｘ）に置き換えて数値積分を行えば、高速な収束が期待できる。なお、このような置き換えを繰り返し行えば、上記の係数の総和がより０に近くなるため、より高速な収束が期待できる関数に被積分関数Ａ（ｘ）を変形していくことができるようになる。
【００９９】
発明者らによるこの研究の成果に基づいて、以下に、本発明の実施形態について説明する。
【０１００】
（第１の実施形態）
図１は、本発明の第１の実施形態に係る数値積分装置の概要構成を示す模式図である。図２は、当該数値積分装置において実行される数値積分方法の処理の流れ示すフローチャートである。以下、これらの図を参照して説明する。
【０１０１】
数値積分装置１０１は、関数Ａ（・）の領域Ｒにおける積分を求め、生成部１０２と、Ｃ計算部１０３と、Δ計算部１０４と、Ｂ計算部１０５と、Ａ更新部１０６と、積分部１０７と、を備える。
【０１０２】
すなわち、生成部１０２は、所定の初期値ｘｌから、漸化式
ｘｊ＋ｌ＝Ｆ（ｘｊ）  （ｊ≧１）
により系列
ｘ１，ｘ２，…
を生成する。
【０１０３】
この関数Ｆとしては、たとえば、チェビシェフ多項式を利用することができる。図３は、チェビシェフ多項式を表したグラフである。チェビシェフ多項式は、整数ａを次数とするとき、
Ｔ（ａ，ｃｏｓθ）＝ｃｏｓ（ａθ）
のように、余弦関数の加法定理により定義することができる。一方、以下のように、有理多項式で直接表現することもできる。
Ｔ（０，ｘ）＝１
Ｔ（１，ｘ）＝ｘ
Ｔ（２，ｘ）＝２ｘ２−１
Ｔ（３，ｘ）＝４ｘ３−３ｘ
【０１０４】
チェビシェフ多項式ｙ＝Ｔ（ａ，ｘ）は、いずれも、開区間−１＜ｘ＜１を開区間−１＜ｙ＜１に写像する有理写像である。
【０１０５】
本図には、次数２から５のチェビシェフ多項式が、ｙ＝Ｔ（２，ｘ），ｙ＝Ｔ（３，ｘ），ｙ＝Ｔ（４，ｘ），ｙ＝Ｔ（５，ｘ）の形式でグラフ表示されている。横軸がｘ軸、縦軸がｙ軸である。
【０１０６】
上述するように、このようにしてＦとして２次以上のチェビシェフ多項式を利用すると、生成される乱数系列の極限分布は、
ρ（ｘ）＝１／〔π（１−ｘ２）１／２〕
となる。
【０１０７】
したがって、まず、生成部１０２は、上記漸化式を使って、Ｍ個の乱数を生成する（ステップＳ２０１）。ここで、後述するように、ステップＳ２０１は繰り返し処理の中の１つのステップであるが、今回の繰り返しで生成された乱数を
ｘ１，ｘ２，…ｘＭ
とおくこととする。なお、次回の繰り返しでは、新しい乱数系列を生成してもよいし、今回と同じ乱数系列を利用してもよい。
【０１０８】
ついで、Ｃ計算部１０３は、生成されたＭ個の乱数を用いて、ラベルｐ∈Ｐ，ｑ∈Ｑ（ｐ）のそれぞれについて、
Ｃｐ，ｑ＝（１／Ｍ）Σｊ＝ｌＭＡ（ｘｊ）φｐ，ｑ（ｘｊ）／ρ（ｘｉ）
を計算する（ステップＳ２０２）。
【０１０９】
そして、計算されたＣｐ，ｑが終了条件を満たすか否かを判断する（ステップＳ２０３）。
【０１１０】
満たさない場合（ステップＳ２０３；Ｎｏ）、Δ計算部１０４は、
Δｐ＝−Σｑ∈Ｑ（ｐ）Ｃｐ，ｑ
を計算する（ステップＳ２０４）。
【０１１１】
ついで、Ｂ計算部１０５は、所定のｑ（Ｐ）∈Ｑ（Ｐ）を用いて、関数
Ｂ（ｘ）＝Ａ（ｘ）＋Σｐ∈ＰΔｐφｐ，ｑ（ｐ）（ｘ）ρ（ｘ）
を計算する（ステップＳ２０５）。ここで得られる関数Ｂ（ｘ）は、当該領域Ｒにおける数値積分としてはＡ（ｘ）と等価な結果が得られるものであり、Ａ（ｘ）よりも収束が速いことが期待されるものである。
【０１１２】
さらに、Ａ更新部１０６は、計算された関数Ｂ（ｘ）で関数Ａ（ｘ）を置き換える（ステップＳ２０６）。
【０１１３】
コンピュータでこれを実現した場合、関数Ａ（ｘ）や関数Ｂ（ｘ）の定義は、コンピュータ上のＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）に格納することとなるため、関数の置き換えはＲＡＭメモリに記憶された情報のコピーによって行うことができる。
【０１１４】
そして、ステップＳ２０１に戻る。
【０１１５】
一方、計算されたＣｐ，ｑが終了条件を満たす場合（ステップＳ２０３；Ｙｅｓ）、生成部１０２は、所定の整数Ｎに対してＮ個の乱数を取得する（ステップＳ２０７）。なお、当該乱数は、ステップＳ２０１で生成されたものと同じ乱数を利用してもよいし、新たな乱数を生成してもよい。そこで、理解を容易にするため、ステップＳ２０７において取得された乱数を
ｘ１，ｘ２，…ｘＮ
とおくこととする。
【０１１６】
そして、積分部１０７は、
Ｉ＝（１／Ｎ）Σｉ＝ｌＮＡ（ｘｉ）／ρ（ｘｉ）
を計算し、Ｉを数値積分結果として（ステップＳ２０８）、本処理を終了する。
【０１１７】
なお、終了条件としては、以下のようなものを採用することができる。
（ａ）計算されたΔｐについて、所定の誤差εに対して
Σｐ∈ＰΔｐ２≦ε
が成立する。
（ｂ）任意のラベルｐ∈Ｐについて、所定の誤差εに対して、
｜Δｐ｜≦ε
が成立する。
【０１１８】
これらの場合は、ステップＳ２０３の判断に先立ってステップＳ２０４の処理を実行する。すなわち、Δ計算部１０４がまずΔｐを計算し、この値を利用して、ステップＳ２０３で終了条件を満たすか否かを判断すればよい。
【０１１９】
領域ＲがＤ次元の領域であるときには、ラベルを以下のように選択することができる。
【０１２０】
すなわち、所定のＤ個の正整数Ｐｍａｘｌ，…，ＰｍａｘＤについて
Ｐ＝｛（ｐ［１］，…ｐ［Ｄ］）｜
ｐ［１］，…，ｐ［Ｄ］は整数；
１≦ｐ［１］≦Ｐｍａｘｌ；…；
１≦ｐ［Ｄ］≦ＰｍａｘＤ｝
とする。
【０１２１】
また、任意のｐ∈Ｐと、所定のＤ個の正整数Ｑｍａｘｌ（ｐ），…，ＱｍａｘＤ（ｐ）と、について
Ｑ（ｐ）＝｛（ｑ［１］，…，ｑ［Ｄ］）｜
ｑ［１］，…，ｑ［Ｄ］は整数；
１≦ｑ［１］≦Ｑｍａｘｌ（ｐ）；…；
１≦ｑ［Ｄ］≦ＱｍａｘＤ（ｐ）｝
とする。
【０１２２】
さらに、力学系を定義する関数Ｆ（・）として、チェビシェフ多項式を採用して、領域Ｒの次元Ｄと、いずれかの整数ａｌ，…，ａＤ（ａｌ≧２，…，ａＤ≧２）と、ｙ＝（ｙｌ，…，ｙＤ）∈Ｒとについて、
Ｆ（ｙ）＝（Ｔ（ａｌ，ｙｌ），…，Ｔ（ａＤ，ｙＤ））
のように定義する。
【０１２３】
この場合、極限密度分布関数ρ（・）は、
ρ（ｘ）＝１／〔πＤ（１−ｙｌ２）１／２…（１−ｙＤ２）１／２〕
となり、
ｐ＝（ｐ［１］，…，ｐ［Ｄ］）；
ｑ＝（ｑ［１］，…，ｑ［Ｄ］）
に対して、基底関数φｐ，ｑ（・）は、
φｐ，ｑ（ｙ）＝２Ｄ／２Ｔ（ｐ［１］ａｌｑ［１］，ｙｌ）…Ｔ（ｐ［Ｄ］ａＤｑ［Ｄ］，ｙＤ）
となる。
【０１２４】
これらの諸元を用いて、数値積分を行うことができる。以下では、数値積分の実験結果を示す。生成部１０２で用いる力学系として
Ｆ（ｘ）＝Ｔ（２，ｘ）
を採用し、被積分関数として
Ａ（ｘ）＝ｘ＋２Ｔ（２，ｘ）＋３Ｔ（４，ｘ）
を採用して、本手法により、数値積分を行い、その結果の誤差の分散を調べた。
【０１２５】
ステップＳ２０１〜Ｓ２０６を適当な回数繰り返した後、ステップＳ２０８のモンテカルロ法において使った乱数の個数Ｎと、数値積分結果の誤算の分散Ｖ（ε）との関係を見ると、
ｌｏｇ（Ｎ）＝２．００⇒ｌｏｇ（Ｖ（ε））＝−３．１５；
ｌｏｇ（Ｎ）＝３．００⇒ｌｏｇ（Ｖ（ε））＝−５．１０；
ｌｏｇ（Ｎ）＝３．５０⇒ｌｏｇ（Ｖ（ε））＝−６．０７；
ｌｏｇ（Ｎ）＝４．００⇒ｌｏｇ（Ｖ（ε））＝−６．８５
となっており、
Ｖ（ε）∝１／Ｎ１．７
が成立しており、通常のモンテカルロ法の収束
Ｖ（ε）∝１／Ｎ
に比べて、格段に高速に収束していることがわかる。
【０１２６】
（その他の実施形態）
本実施形態は、上記の数値積分装置１０１を複数用意して並列に実行させるものである。図４は、本実施形態の概要構成を示す模式図である。
【０１２７】
本実施形態の数値積分装置４０１は、初期値付与部４０２と、平均出力部４０３と、複数の上記数値積分装置１０１とを備える。
【０１２８】
すなわち、初期値付与部４０２は、各数値積分装置１０１の生成部１０２に、互いに異なる初期値を当該所定の初期値ｘｌとして与える。すると、これらの数値積分装置１０１は、並列に数値積分を実行して、それぞれの積分部１０７から数値積分の結果Ｉを出力する。
【０１２９】
本実施形態では、互いに異なる初期値のそれぞれに対して、各数値積分装置１０１の生成部１０２、Ｃ計算部１０３、Δ計算部１０４、Ｂ計算部１０５、Ａ更新部１０６、ならびに、積分部１０７により行われる処理が、並列に行われることとなる。
【０１３０】
そして、平均出力部４０３は、それぞれの数値積分装置１０１が出力した結果の相加平均を、数値積分結果として出力する。
【０１３１】
本実施形態の場合、終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算部１０４により計算されるΔｐの自乗和の総和が、所定の誤差ε以下である場合に満たされるようにすることができるほか、終了条件は、互いに異なる初期値のそれぞれを与えた結果Δ計算部１０４により計算されるΔｐのいずれもが、所定の誤差ε以下である場合に満たされるようにすることができる。
【０１３２】
これにより、より精度の高い数値積分を高速に行うことができると期待される。
【０１３３】
なお、複数の数値積分装置１０１を用意して並列に動作させるのではなく、１つ用意して逐次に実行させることによっても、同様の処理を行うことが可能である。
【０１３４】
【発明の効果】
以上説明したように、本発明によれば、収束の高速な数値積分装置、数値積分方法、ならびに、これらをコンピュータ上で実現するプログラムを提供することができる。
【図面の簡単な説明】
【図１】本発明の数値積分装置の実施形態の概要構成を示す模式図である。
【図２】数値積分装置において実行される処理の流れを示すフローチャートである。
【図３】チェビシェフ多項式の様子を示すグラフである。
【図４】本発明の他の実施形態の概要構成を示す模式図である。
【符号の説明】
１０１  数値積分装置
１０２  生成部
１０３  Ｃ計算部
１０４  Δ計算部
１０５  Ｂ計算部
１０６  Ａ計算部
１０７  積分部
４０１  数値積分装置
４０２  初期値付与部
４０３  平均出力部
【図１】
0003711270.tif 000002
【図２】
0003711270.tif 000003
【図３】
0003711270.tif 000004
【図４】
0003711270.tif 000005