広義積分とは
無限大の積分である。数学的には
\int_{-\infty}^{\infty}f(x)dx=\lim_{a\rightarrow\infty}\int_{-a}^{a}f(x)dx
と極限で定義される。数値積分もそれに乗っ取り、有限範囲からで積分したものを範囲を広げつつ、その差が誤差以下になれば打ち切りる。
しかし、普通の数値積分を使うと積分範囲が大きくなりすぎ数値的な丸めや誤差が大きくなり収束しない。そこで次のような変数変換を使う。
2重指数関数型変換積分
前述の発散を抑えるため変数変換で影響のある範囲を0付近に集める。この変数変換に双曲関数$sinh,consh,tanh$を用いたものが2重指数関数型変換積分である。双曲関数が指数型で書かれているのでこの名称がつている。$1/x$のような発散を含む積分も変数変換で取り除くので安定しやすい。有限範囲でも変換後は台形で近似しやすい形になるので使いやすさが上がる(らしい)。
| 関数 | 変換前範囲 | 変換後範囲 | 使用する例 |
|---|---|---|---|
| tanh((π/2)sinh(t)) | -1:1 | -∞:∞ | 有限積分,２つの極の間の積分 |
| sinh((π/2)sinh(t)) | -∞:∞ | -∞:∞ | 無限領域の積分 |
| exp((π/2)sinh(t)) | 0:∞ | -∞:∞ | 半無限領域の積分 |
| exp(t-exp(t)) | 0:∞ | -∞:∞ | f(x)exp(-x)で書ける時の特例(?) |
要は変数変換なのでこう書ける。
#pragma once
#include "Integrator.hh"
#include "Differentaitor.hh"
namespace MyMath{
namespace ImproperIntegrator{
double DE_mInf_Inf(const std::function<double(double)> &f, int N=10000, double error=DEFUALT_ERROR);
double DE_0_Inf(const std::function<double(double)> &f, int N=10000, double error=DEFUALT_ERROR);
};
};
$x=f(t)$ならば$dx=f'(t)dt$なので非積分関数にconverterの微分がかかる。
変換関数は$~7$付近で急速に値が大きくなる(倍精度でもオーバーフローする)。のでその半分からスタートして0.1幅で刻んでいる。
#include "ImproperIntegrator.hh"
using namespace MyMath;
double ImproperIntegrator::DE_mInf_Inf(const std::function<double(double)> &f, int N, double error){
double width=3.5;
auto converter=[](double x){ return sinh(PI*sinh(x)/2.0); };
auto integral_func=[f, converter](double x){ return f(converter(x))*Differentaitor::central(converter, x); };
double val=Integrator::trapezoidal(integral_func, -width, width, N);
while( !std::isnan(val) && !std::isnan(width) ){
width+=0.1;
double next=Integrator::trapezoidal(integral_func, -width, width, N);
if( fabs(next-val)<error ) return next;
val=next;
}
throw std::runtime_error("ImproperIntegrator::DE_mInf_Inf can not calculate");
}
double ImproperIntegrator::DE_0_Inf(const std::function<double(double)> &f, int N, double error){
double width=3.5;
auto converter=[](double x){ return exp(PI*sinh(x)/2.0); };
auto integral_func=[f, converter](double x){ return f(converter(x))*Differentaitor::central(converter, x); };
double val=Integrator::trapezoidal(integral_func, 0, width, N);
while( !std::isnan(val) && !std::isnan(width) ){
width+=0.1;
double next=Integrator::trapezoidal(integral_func, 0, width, N);
if( fabs(next-val)<error ) return next;
// std::cout<<width<<" "<<val<<" "<<next<<std::endl;
val=next;
}
throw std::runtime_error("ImproperIntegrator::DE_0_Inf can not calculate");
}
広義積分からπが計算できた
cout<<"S1/(x*x+1)dx [-inf:inf] = "<<ImproperIntegrator::DE_mInf_Inf([](double x){ return 1/(x*x+1); })<<endl;
return 0;
}
$\frac{1}{x^2+1}$の積分で$\pi$になる(留数定理で計算できます)。
その他$1/x$のような発散を含むものも発散点近くの影響をちらして積分するので精度があがるらしい。一般の数値計算の教科書でもここまで書いているのは珍しい気がする。
その他の変数変換も稀に用いられることがあるらしいが原理は同じなのでconverterを変えるだけで対応できるだろう。
ただし対応限界はあり$sin,cos$などの振動しているものは収束しません。
参考にしたサイト
Wikipedia:二重指数関数型数値積分公式
数値積分（7）（DE公式 続き）
数式処理ソフトで遊ぼう（DERIVE デライブ）へようこそ---上の上位ページ