teratail header banner
teratail header banner
質問するログイン新規登録

回答編集履歴

1

暫定回答追加

2016/09/24 02:50

投稿

cesolution
cesolution

スコア217

answer CHANGED
@@ -1,2 +1,55 @@
1
1
  sukiyaki様
2
- 一通りコードを確認しましたが、少なくともシンプソン積分、2分法のコーディングは問題無さそうです。ちょっと被積分関数の形状がパッとは分からないのですが、恐らくかなり複雑な形状をしており、0~1の間に複数解を持つのではないでしょうか?今回は答えが分かっているようなので、2分法で挟む範囲を狭めてみてはどうでしょうか?
2
+ 一通りコードを確認しましたが、少なくともシンプソン積分、2分法のコーディングは問題無さそうです。ちょっと被積分関数の形状がパッとは分からないのですが、恐らくかなり複雑な形状をしており、0~1の間に複数解を持つのではないでしょうか?今回は答えが分かっているようなので、2分法で挟む範囲を狭めてみてはどうでしょうか?
3
+
4
+ 2016/9/24追記
5
+  いくつか試してみましたが、現状の関数のままだと、関数の出力が初期値を入力した段階で#INDになってしまうようです。
6
+ 試しに被積分関数をtanh(x/(2*T))にしてみても同様でした。
7
+  ひとまず収束して解を得られるコードを以下、一例としてご紹介します。(スマートなコードではないですが、、、)
8
+ mainの中で以下のbisection2()を実行していただければと思います。
9
+
10
+ ```C++
11
+
12
+ double integrand2(double x, double T){
13
+ return tanh(x/(2*T))/max(0.000001,x);
14
+ };
15
+
16
+ double int_simpson2(int N, double a, double b, double T){
17
+ double dx=(b-a)/N;
18
+ double sum1=0.,sum2=0.;
19
+ for (int i=1; i<N/2; i++) {
20
+ sum1+=integrand2(a+2*i*dx, T);
21
+ sum2+=integrand2(a+(2*i-1)*dx, T);
22
+ }
23
+ sum2+=integrand2(a+(N-1)*dx, T);
24
+ return dx*(integrand2(a, T)+integrand2(b, T)+2*sum1+4*sum2)/3.;
25
+ }
26
+
27
+ void bisection2(){
28
+ int count=0;
29
+ double l=0.;
30
+ //double r=1.0;
31
+ double r=0.15;
32
+ double VN=0.4;//相互作用
33
+ double m;
34
+ do{
35
+ count++;
36
+ m=(l+r)/2.0;
37
+ if ((int_simpson2(100,0.,1.,m)-1.0/VN)*(int_simpson2(100,0.,1.,l)-1.0/VN)<0) r=m;
38
+ else l=m;
39
+ if(count==20){
40
+ cout << "収束しませんでした";
41
+ exit(1);
42
+ }
43
+ //cout << int_simpson2(100,0.,1.,m)-1.0/VN << endl;
44
+ //cout << m << endl;
45
+ }while (!(fabs(int_simpson2(100,0.,1.,m)-1.0/VN)<0.001));
46
+ cout << m << endl;
47
+ }
48
+
49
+ ```
50
+
51
+ なお、2分法の範囲を狭めていますが、理由は以下になります。
52
+
53
+ ![イメージ説明](8969fd648079968432c9badb9ac3fc1b.png)
54
+
55
+ ![イメージ説明](de5452d0acd5976e83ab5cc50ce610d7.png)