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

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

新規登録して質問してみよう
ただいま回答率
85.50%
C++

C++はC言語をもとにしてつくられた最もよく使われるマルチパラダイムプログラミング言語の1つです。オブジェクト指向、ジェネリック、命令型など広く対応しており、多目的に使用されています。

Q&A

解決済

2回答

2794閲覧

積分と二分法

sukiyaki

総合スコア15

C++

C++はC言語をもとにしてつくられた最もよく使われるマルチパラダイムプログラミング言語の1つです。オブジェクト指向、ジェネリック、命令型など広く対応しており、多目的に使用されています。

0グッド

1クリップ

投稿2016/09/22 00:24

###前提・実現したいこと
ギャップ方程式での臨界温度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
idx);
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/ツール等のバージョンなど)

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答2

0

ベストアンサー

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

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

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

イメージ説明

イメージ説明

投稿2016/09/22 02:51

編集2016/09/24 02:50
cesolution

総合スコア217

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

sukiyaki

2016/09/23 01:05

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

2016/09/24 09:07

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

0

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

double integrand(double x){

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

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

投稿2016/09/22 00:40

編集2016/09/22 00:41
strike1217

総合スコア651

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

sukiyaki

2016/09/22 01:40

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問