質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.83%

積分と二分法

解決済

回答 2

投稿

  • 評価
  • クリップ 1
  • VIEW 1,759

sukiyaki

score 15

前提・実現したいこと

ギャップ方程式での臨界温度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+2*i*dx);
sum2+=integrand(a+(2*i-1)*dx);
}
sum2+=integrand(a+(N-1)*dx);
return dx*(integrand(a)+integrand(b)+2*sum1+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/ツール等のバージョンなど)

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 2

checkベストアンサー

+1

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

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

double integrand2(double x, double T){
return tanh(x/(2*T))/max(0.000001,x);
};

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

void bisection2(){ 
int count=0; 
double l=0.; 
//double r=1.0;
double r=0.15; 
double VN=0.4;//相互作用
double m; 
do{ 
count++; 
m=(l+r)/2.0; 
if ((int_simpson2(100,0.,1.,m)-1.0/VN)*(int_simpson2(100,0.,1.,l)-1.0/VN)<0) r=m; 
else l=m; 
if(count==20){ 
cout << "収束しませんでした"; 
exit(1); 
} 
//cout << int_simpson2(100,0.,1.,m)-1.0/VN << endl; 
//cout << m << endl; 
}while (!(fabs(int_simpson2(100,0.,1.,m)-1.0/VN)<0.001)); 
cout << m << endl; 
}

なお、2分法の範囲を狭めていますが、理由は以下になります。

イメージ説明

イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2016/09/23 10:05

    二分法の範囲をかなり狭めて試してみましたが、答えが最初の一桁しか合いませんでした。
    試してみたことにあるように被積分関数をTにしてみて計算した結果も合わなかったので、二分法の関数fがうまく渡せてないのかなと思うのですが...

    キャンセル

  • 2016/09/24 18:07

    非常にわかりやすい回答ありがとうございます。おかげで期待していた答えに近い値がでました。まだ経験が乏しくなかなか原因を見つけるのが難しかったので助かりました。

    キャンセル

  • 2016/09/25 17:34

    天才ですね!!

    キャンセル

0

間違っていましたら、ゴメンナサイ。

double integrand(double x){

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

の部分なんですが、式から推測するに "tanh(x/(2*T))/x-1"にしてみてはどうでしょ?
integrand(x-1);にしても数式と一致しません。

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2016/09/22 10:40

    式がわかりにくくてすいませんが、結局求めたいのは
    積分=1/(VN)
    の解なので被積分関数はあっていると思います。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.83%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る