積分と二分法

###前提・実現したいこと
ギャップ方程式での臨界温度Tを二分法とシンプソン積分を使って求めたい。
解きたい方程式
int^{1}_{0}dx tanh(x/(2*T))/x-1/(VN)=0
###発生している問題・エラーメッセージ

T=0.09程度の値が期待されるが、計算結果は0.9程度になってしまった

###該当のソースコード
C++

#include <iostream>
#include <iomanip>
#include <cmath>
#include <complex>
#define eps 1.0e-5
using namespace std;

double T;
double integrand(double x){

return tanh(x/(2*T))/x;

};

double int_simpson(double (f)(double), int N, double a, double b){
double dx=(b-a)/N;
double sum1=0.,sum2=0.;
for (int i=1; i<N/2; i++) {
sum1+=integrand(a+2idx);
sum2+=integrand(a+(2i-1)dx);
}
sum2+=integrand(a+(N-1)dx);
return dx(integrand(a)+integrand(b)+2sum1+4*sum2)/3.;
}

double f(double){
double VN=0.4;//相互作用
return int_simpson(integrand,100,0.,1.)-1./(VN);
}

void bisection(){
int count=0;
double l=0.;
double r=1.;
double m;
do{
count++;
m=(l+r)/2.0;
if (f(m)*f(l)<0) r=m;
else l=m;
if(count==100){
cout << "収束しませんでした";
exit(1);
}
}while (!(fabs(l-r)<eps));
cout << m << endl;
}

int main (){

bisection();

return 0;

}

###試したこと
被積分関数を簡単にTとして計算した場合も期待される答えが得られなかった。（二分法の区間を[1,3]とした）

###補足情報(言語/FW/ツール等のバージョンなど)

行動規範の内容に同意します

回答2件

ベストアンサー

sukiyaki様
一通りコードを確認しましたが、少なくともシンプソン積分、２分法のコーディングは問題無さそうです。ちょっと被積分関数の形状がパッとは分からないのですが、恐らくかなり複雑な形状をしており、0～1の間に複数解を持つのではないでしょうか？今回は答えが分かっているようなので、２分法で挟む範囲を狭めてみてはどうでしょうか？

2016/9/24追記
いくつか試してみましたが、現状の関数のままだと、関数の出力が初期値を入力した段階で#INDになってしまうようです。
試しに被積分関数をtanh(x/(2*T))にしてみても同様でした。
ひとまず収束して解を得られるコードを以下、一例としてご紹介します。（スマートなコードではないですが、、、）
mainの中で以下のbisection2()を実行していただければと思います。

C++
1
2double integrand2(double x, double T){
3return tanh(x/(2*T))/max(0.000001,x);
4};
5
6double int_simpson2(int N, double a, double b, double T){ 
7double dx=(b-a)/N; 
8double sum1=0.,sum2=0.; 
9for (int i=1; i<N/2; i++) { 
10sum1+=integrand2(a+2*i*dx, T); 
11sum2+=integrand2(a+(2*i-1)*dx, T); 
12} 
13sum2+=integrand2(a+(N-1)*dx, T); 
14return dx*(integrand2(a, T)+integrand2(b, T)+2*sum1+4*sum2)/3.; 
15}
16
17void bisection2(){ 
18int count=0; 
19double l=0.; 
20//double r=1.0;
21double r=0.15; 
22double VN=0.4;//相互作用
23double m; 
24do{ 
25count++; 
26m=(l+r)/2.0; 
27if ((int_simpson2(100,0.,1.,m)-1.0/VN)*(int_simpson2(100,0.,1.,l)-1.0/VN)<0) r=m; 
28else l=m; 
29if(count==20){ 
30cout << "収束しませんでした"; 
31exit(1); 
32} 
33//cout << int_simpson2(100,0.,1.,m)-1.0/VN << endl; 
34//cout << m << endl; 
35}while (!(fabs(int_simpson2(100,0.,1.,m)-1.0/VN)<0.001)); 
36cout << m << endl; 
37}
38