暇つぶしで複素積分の近似計算に挑戦

---

仕事に飽きたので（最近よく飽きるのだ）眠気覚ましと頭の体操の為に適当な数学の問題に挑戦する。
ノリはこれに近い。
挑戦する問題はこれ↓だ！（Windows標準の数式エディタで書いたので細部がいまいちだがまあいいや）

---

---

 

この問題を選んだことに意味はない。
ふと頭の中に浮かんだ数式だったというだけである。
強いて言うなら冒頭紹介した「ノリはこれに近い」の問題でSin(x)を扱ったというのから出てきたところはあるかもしれない。

被積分関数は、x→∞の極限で0に収束するだろうことは容易に想像できる。
よってこの積分値自身もどこかに収束すると予想できる。
なお、この問題、理論的に解く場合、複素積分の考え方が必要になるらしい。
ググったらすぐ出てくるほど複素積分の問題としてはオーソドックスなもののようだ。
⇒たとえばこことか

そういえば大学の時「物理数学」（物理のための数学テクや知識を学ぶ）とかいう講座で習った気がするが
最早何書いてあるかわからんほど意味不明な展開なので
理論的なアプローチは完全に諦めてプログラミングによるアプローチに切り替えていく。
（まあ最初からそのつもりだったわけだが）
ちなみに、上記の解説によれば、この問題の解はだそうだ。
果たして「プログラミングによるアプローチ」でこれに近い値に迫れるだろうか？

---

相変わらずやり方は昔と同じで、「細かく刻んで足しこむ」というリーマン和を使う。
刻み幅×関数値で面積を計算し、それを延々と足しこんで合計値を求め、疑似的に積分を実現させる。
「∞」が終端なのでできればwhile(true)でも使って延々と計算を続けたいところだが、
（やってみればわかるが）分母のパワーが強すぎるため非積分関数の収束は結構早く、
例えば1000まで進むと最早0付近をうろちょろしだして、変化という変化がなくなる。
このため積分値もこの辺で大体収束に近づいてくるといえる
（もちろん延々と計算を続ければもっと精度の高い積分値に迫ることは可能なのだろうが）
なのでとりあえず100,000あたりを狙っていこうと思う。

ちなみに今回、分母にxがいる関係で、x=0は初期値として取れない（いわゆる#DIV/0!が起きる）
このため「0に近い値」を取ることにする。
刻み幅を含め、とりあえず0.0001あたりを採用してみる。

この問題のポイントは、解がになるという点だ。
つまり円周率の半分なのだ。
逆に言えば計算して求めた合計値を×２すると円周率の近似値になっているはずだ。
「理論解に近づいた」ことを確認するよりは
「２倍した値が円周率に近づいた」ことを確認する方がぱっと見わかりやすいので、
プログラムの計算過程では合計値×２も計算して標準出力させていくことにする。

BigDecimalでも使えばもっと精度のいい数値が出るのだろうが、面倒なのでdoubleだけで実装する。

import java.math.*;  
import java.util.*;  
import java.text.*;  
public class Himatsubushi {
private static final double ROCK_END_ROLL = 99999.9;  
  
private static final DateFormat LOG_FORMAT = new SimpleDateFormat("yyyy/MM/dd HH:mm:ss.SSS");  
private static final int REST = 10000000;  
  
public static void main(String[] args) {  

	exec();  
}  
  
private static void exec() {  
	printLog("★開始");  
	double ans = 0.0;  
	double pi = 0.0;  
	double n = 0.0;  
	boolean isContinue = true;  
	int round = 0;  
	  
	double sabun = 0.0001; // 初期値にも使う  
	do {  
		double nWk = n == 0.0 ? sabun : n;  
		double calc = Math.sin(nWk) / nWk;  
		ans = ans + (calc * sabun);  
		pi = ans * 2.0;  
		  
		round++;  
		n = n + sabun;  
		if (round % REST == 0 && round > 0) {  
			printLog(round + "回目終了...合計値＝" + ans + "、合計値×２（π想定）＝" + pi + "、現在値＝" + calc);  
		}  
		isContinue = n <= ROCK_END_ROLL;  
		  
	} while(isContinue);  
	pi = ans * 2.0;  
	printLog("＊最終合計値＝" + ans + "、最終合計値×２（π想定）＝" + pi);  
	  
	printLog("★終了");  
}  
  
  
private static void printLog(Object msg) {  
	StringBuilder sb = new StringBuilder();  
	sb.append(LOG_FORMAT.format(new Date()));  
	sb.append(" ");  
	if (msg != null) {  
		sb.append(msg.toString());  
	}  
	System.out.println(sb.toString());  
  
}  

}

ちなみにこの結果、10億回ループして計算して、最終的には

No種類計算値理論値

1	合計値（π/2想定）	1.5708517123342043	1.5707963267949
2	合計値×２（π想定）	3.1417034246684086	3.14159265358979

となった。

う～ん、なんかちょっとおしいね！
欲を言えばもうちょっとお近づきになりたかったところだが。
まあ暇つぶしでやるにはこんなもんかあ。

ちなみに、Sin(x)は-1～1の間を行ったり来たりする関数なので、
非積分関数の値は正の数になったり負の数になったりする。（まあxがでかくなると正負に関わらず0付近にいるが）
このためタイミングによっては↑の結果に比べて理論値に近い値を取ることも有る。
例えば1億2千万回目くらいの計算結果では

No種類計算値理論値

1	合計値（π/2想定）	1.5707934936780714	1.5707963267949
2	合計値×２（π想定）	3.141586987356143	3.14159265358979

という結果を得ている。
↑の結果に比べると理論値に近いことがわかる。

---