よくある質問というよりは印象にのこている質問です。
- Q0 : 積分 integral 1/x dx from 0 to 1 が計算できない。 (作り話のように思うかも知れませんが、 本当に聞いてくる人がいます。)
- これは発散する積分です。計算できません。 しかし積分
integral sin(x) dx from 0 to infinity
のような収束しない振動積分は計算ができて (intdeo使用のときのみ)積分
integral sin(x)*exp(-eps*x) dx from 0 to infinity
の eps が 0 の極限の値を返します。
- Q1 : 積分 integral 1/sqrt(1-x^2) dx from 0 to 1 が計算できない。
- 被積分関数が発散する点 x = 1 近傍で桁落ち(正確には 積分パッケージ内での情報落ち + 被積分関数での桁落ち) がおきています。要するに計算機内部で
1/sqrt(1-0.999999999999999999999xxx) = 1/0
としてしまいます。 被積分関数が発散している x = 1 近傍での被積分関数の 積分値への寄与は非常に大きいため Not a Number や誤差の大きい値を返されると計算不可能になります。 いくつかの解決策の例を以下に示します。
- 被積分関数が発散する点を原点にずらし
integral 1/sqrt(-x*(2+x)) dx from -1 to 0
と桁落ちがおきないように変形する。
integral 1/sqrt(1-(1+x)^2) dx from -1 to 0
のままでは意味がないので注意!- 区間を変えたくない場合、あるいは積分の両端が 特異点の場合はパッケージを修正して二つの引数 x および 1-x (あるいは必要な方)を被積分関数に 渡すようにする。 さらに被積分関数は桁落ちがおきないように 修正しなければならない。
- 低い精度(たとえば４桁)でいいという場合は
integral 1/sqrt(1-x^2) dx from 0 to 0.999999999
のように必要な精度内で区間を狭めれば変形無しで 直接計算できる。- Q2 : 積分 integral sqrt(x)/(exp(x)-1) dx from 0 to 1 が計算できない。
- 原因はQ1とほぼ同じです。x = 0 で桁落ちがおきています。
対策は exp(x)-1 = 2*exp(x/2)*sinh(x/2) と変形します。
exp(x)-1 を桁落ち無しで返す関数が使える場合は 何も問題ありません。
- Q3 : 積分 integral 1/sqrt(log(x)) dx from 1 to 2 が 計算できない。
- 原因はQ1とほぼ同じです。x = 1 で桁落ちがおきています。
対策は積分区間をずらして
integral 1/sqrt(log(1+x)) dx from 0 to 1
として log(1+x) = 2*atanh(x/(2+x)) と変形します。 atanh(x) が使えない場合は x = 0 近傍での log(1+x) のテイラー展開より計算します。
そんな面倒なことはしたくないという場合は変数変換 x = exp(t) を行い
integral exp(t)/sqrt(t) dt from 0 to log(2)
により計算します。
log(1+x) を桁落ち無しで返す関数が使える場合は 何も問題ありません。
- Q4 : 積分 integral log(x)/(x-1) dx from 0 to 1 が計算できない。
- 原因はQ1とほぼ同じですがこの場合は x = 1 が除去可能な 特異点で発散がないので簡単に解決できます。
- 対策その１
積分区間を 0 から 1-eps にする。 eps は許容相対誤差である。- 対策その２
被積分関数を x = 1 での極限値 1 を返すようにする。- Q5 : 積分 integral 1/cosh(x^2) dx from 0 to infinity が 計算できない。
- 原因は cosh(x^2) がオーバーフローすることです。
intdei では最初に x = 1 を中心に、許容相対誤差の 平方根からその逆数までの区間での被積分関数の挙動を調べます。 その時点でこのオーバーフローが起きます。
IEEE754 double 準拠の場合、オーバーフロー限界は 1.797e308 で exp(x), cosh(x), sinh(x) の引数の限界は 709 です。 この場合、無限大に相当する値がサポートされていれば オーバーフローが起きても正しく計算できる場合があります。 しかし調べた限りではオーバーフローやアンダーフローの 扱いは機種やコンパイラ、コンパイラの最適化オプションによって すべて異なるようです。
いくつかの解決策の例を以下に示します。
- 関数がオーバーフローする引数の値を調べて その値以上で計算せずに 0 を返すようにする。
- 有限区間の積分パッケージ(intde)で計算する。
integral 1/cosh(x^2) dx from 0 to infinity
は倍精度の場合、区間 (0,7) の有限区間積分と 見なして問題はない。
減衰が速く近似的に有限区間積分と見なせる積分は intdei を使うよりも intde を使うほうが効率がよい。- 被積分関数の減衰を遅くできればオーバーフロー は起きなくなる。
減衰を遅くする変数変換 x = sqrt(log(1+t)) を行い
integral 1/(1+(1+t)^2)/sqrt(log(1+t)) dt from 0 to infinity
と変形し、減衰を代数関数的オーダーにすれば intdei によりオーバーフローせずにしかも効率よく計算できる (log(1+t)の桁落ちに注意)。- 指数関数減衰専用のDE公式を使う。
これは効率のよい方法ですが、あまり必要だとは 思わなかったのでこのパッケージには入れてません。 必要ならば 参考文献 を調べて自分で作成してください。 intdei を少し変形するだけでも簡単に作成できます。- Q6 : 積分 integral sqrt(abs(0.7-x)) dx from 0 to 1 が計算できない。
- 被積分関数が端点を除く積分区間内 x = 0.7 で微分不可能です。この場合 x = 0.7 で積分を二つに分割すれば計算できます。
微分不可能や不連続な点がわからないときには 適応型積分ルーチン(たとえば netlib/quadpack ) を用いる方法があります。
- Q7 : 積分 integral 1/(0.1^10+(0.7-x)^2) dx from 0 to 1 が計算できない。
- 被積分関数が端点を除く積分区間内 x = 0.7 で鋭いピークを持っています。この場合 x = 0.7 で積分を二つに分割すれば計算できます。
この場合 Q6 と同様に適応型積分ルーチンを用いる 方法があります。
- Q8 : 積分 integral sin(x)*log(x) dx from 0 to 100000 が計算できない。
- これは振動積分です。
intde で計算する場合ルーチン内の mmax の値を大きくすれば無理に計算させることもできます。 この積分をまじめに計算する場合、どんな積分公式でも 少なくとも振動の回数以上の標本点数が必要です (直感的に明らか)。
しかしこの積分の計算には裏技があって
integral sin(x)*log(x) dx from 0 to 100000
を無限遠点で区間分割して
integral sin(x)*log(x) dx from 0 to infinity -
integral sin(x)*log(x) dx from 100000 to infinity
として振動積分ルーチン intdeo を使うという方法があります。
- Q9 : 積分 integral sin(1/x) dx from 0 to 1 が計算できない。
- これは振動積分です。変数変換 x = 1/t により
integral sin(t)/t^2 dt from 1 to infinity
と変形して intdeo を使います。
- Q10 : 積分 integral sin(x+1/x)/sqrt(x) dx from 0 to infinity が計算できない。
- これは振動積分です。x = 0 と x = infinity の両方で振動しています。
まず積分区間を二つに分けます。
integral sin(x+1/x)/sqrt(x) dx from 0 to 1 +
integral sin(x+1/x)/sqrt(x) dx from 1 to infinity
次に最初の積分を変数変換 x = 1/t で
integral sin(t+1/t)*sqrt(t)/t^2 dt from 1 to infinity +
integral sin(x+1/x)/sqrt(x) dx from 1 to infinity
と変形して intdeo を使います(振動数 omega はこの場合 1 とします)。
- Q11 : Cauchy の主値積分を計算したい。
- 桁落ちを無視すれば普通の積分パッケージでも 強引に計算させることもできるのですが基本的には 別のパッケージが必要(たとえば netlib/quadpack )です。
- Q12 : 漸近的に三角関数以外の振る舞いをする 振動積分あるいは複雑な周期の振動積分 を計算したい。
- これにはうまい方法はなく個別の対応になります。
- Q13 : このパッケージは自動積分ですか。
- 自動積分です。しかし適応型ではありません。
誤差の指定方法は相対誤差で行います。ただし被積分関数の 和の計算で桁落ちが起きた場合、落ちた桁数分の精度は 除外されます。
推定誤差は絶対誤差を返し、 指定された精度に達しないと負で返します。
パッケージ内のパラメータ hoff は初期計算の刻み幅に比例する値、パラメータ mmax は最大標本点数に比例する値になります。
- Q14 : 積分 integral log(x) dx from 0 to exp(1) の計算を 15桁(許容相対誤差 = 1.0e-15)で計算させたが 真の値 0 に対して 3.2e-16 という一桁も一致しない 答えを返す。
- これは当然です。ルーチンでの許容相対誤差 eps は 本当の相対誤差 :
(絶対誤差) / (積分値)
ではなくて
(絶対誤差) / (絶対値をとった被積分関数の積分値)
を意味します。パッケージをいじれば本当の相対誤差に 変更することもできますがそうするとゼロになる積分で Mathematicaの数値積分ルーチンのように深い計算ループに陥って 答えを返さなくなるので注意してください。
- Q15 : 積分の計算で指定した許容相対誤差 eps よりも返される推定誤差 err が大きくなることがあるのはなぜですか。
- これは仕様です。最初のバージョンで、 log(eps) を大まかな桁数、 err を誤差の目安/エラーメッセージとし、 その後の修正で err を実際の誤差よりもなるべく 下回らないようにしたためです。
- Q16 : 積分の計算を大量に行った場合たまに(数百万から 数億回に一回ぐらいの割合で)推定誤差よりも 悪い結果を返す。
- これは当然です。どんな積分ルーチンでも起こることで 確率的に小さくすることはできても完全にゼロにすることは できません。なぜならばたとえば積分公式の標本点で ゼロになり積分値がゼロにならない多項式が簡単に 作れるからです。
このパッケージでは離散化誤差の判定を二重に行って (積分値 I 以外にDE変換後の被積分関数を cosh(t) で割った誤差推定判定以外に意味のない積分 Ir を計算している) この確率を小さくしています。実験より離散化誤差の推定を 誤った場合の実際の誤差は指定した桁数の半分を 下回ることはないようです。
離散化誤差/打切り誤差の判定ミスに対する対策は プログラム中のパラメータ esf を小さくする (たとえば許容相対誤差の平方根程度)ことで改善できます。 また刻み幅パラメータ hoff を小さくすることで 最初のステップでの被積分関数の読み飛ばしを改善できます。 しかしこれらの変更をすると計算時間が遅くなるので 注意してください。
- Q17 : intdeo では等間隔でない周期で振動する積分が 計算できるのですがどういうアルゴリズムを 用いているのですか。
- 正確にはDE公式に Euler変換の加速ルーチンをつけ足す方法を 用いています。
最初のステップで被積分関数が g(x)*sin(omega*x) という形をしているものと仮定してDE変換を行います。 その後、標本点が被積分関数のゼロ点に近づかない場合に 加速を行うようにしています。
次のステップで変換を固定した状態で刻み幅を半分にします。
この反復を誤差が許容値に達するまで繰り返します。 効率をよくするため、 この反復はたいての場合２回ですむようにしてあります。
このルーチンで用いたDE変換は
phi(t) = (t+sqrt(t^2+exp(p-k*cosh(t))))/2
という変換を用いていて、参考文献で用いている変換
phi(t) = t/(1-exp(-2*pi*sinh(t)))
とは形が異なるので(ソースを読む気があるという場合) 注意してください。
DE公式はかなり強力な積分公式ですが、無限区間の積分や 端点で発散する積分を計算する場合は以上のような注意が必要です (DE公式以外でも同様です)。 一番大きな実数が 1.797e308 だと思っている計算機で 極限の計算をさせるのだから当然です。 もしこのような工夫が面倒だと思う場合は多倍長 (たとえば UBASIC )で計算する ことをすすめます。ほんの少しの変更でこれらの多くの問題を 解決できます。
そのうち大幅にパッケージを書き換えようと思っていますが 時間がないのであまり期待しないでください。