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

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

ただいまの
回答率

88.06%

2分法で等高線を求めたものの、期待通りの答えにならない。

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 1,222

score 89

前提・実現したいこと

物理学が専門の大学院生です

二変数関数f(x,y)=cを満たす(x,y)をたくさん求めるプログラムをC言語を用いて、二分法で考えました。しかし、f(x,y)の等高線を数式処理ソフトで書かせた場合と一部一致しません。

発生している問題・エラーメッセージ

左の赤いラインが二分法によって求めたf(x,y)=1.54をみたす点(x,y)をプロットしたもの。一方で、右の青い図が数式処理ソフトによって出力された関数f(x,y)の等高線。二分法でもとめた点(x,y)の一部がかけているのはなぜでしょうか。

イメージ説明
イメージ説明

該当のソースコード

include <stdio.h>
include <math.h>
double f(double,double);

void main(void){
  double x,xmin,xmax,xmid,y,eps;
  int cnt,cntlimit,i,m;
  FILE *file;

  eps=0.001;
  cntlimit=100000;

  file = fopen("data.dat","w");
  cnt=0;

  for(i=-500;i<=500;i++){

    y=2.0+i/100.0;
    xmax=2.5;
    xmin=0.0;

  do{
    cnt ++;
    xmid=(xmax+xmin)/2.0;

    if( f(xmin,y)*f(xmid,y)<0.0) xmax=xmid ;
    else xmin=xmid;

  }while( fabs(xmax-xmin)>eps && cnt < cntlimit );

   if(cnt==10000) fprintf(file,"It doesnt converge \n");
   else if( fabs(f(xmid,y))>0.01 ) fprintf(file,"0\t0\n");
   else fprintf(file,"%f\t%f\n",xmid,y);


  }

  for(i=-500;i<=500;i++){

    y=2.0+i/100.0;
    xmax=-0.0;
    xmin=-2.5;

  do{
    cnt ++;
    xmid=(xmax+xmin)/2.0;

    if( f(xmin,y)*f(xmid,y)<0.0) xmax=xmid ;
    else xmin=xmid;

  }while( fabs(xmax-xmin)>eps && cnt < cntlimit );

   if(cnt==10000) fprintf(file,"It doesnt converge \n");
   else if( fabs(f(xmid,y))>0.01 ) fprintf(file,"0\t0\n");
   else fprintf(file,"%f\t%f\n",xmid,y);


  }


  fclose(file);
}

double f(double kx,double ky){
  double c=1.98;
  double w2p,w2d,eE,t,s,d,Ec,Ev;

  t=-1;
  s=0;
  d=1;
  eE=0.5;

  w2p=3+4*cos(sqrt(3)*(kx+eE)/2.0)*cos(ky/2.0)+2*cos(ky);
  w2d=3+4*cos(sqrt(3)*(kx-eE)/2.0)*cos(ky/2.0)+2*cos(ky);

  Ec=(2*t*s*w2p+sqrt(d*d-d*d*s*s*w2p+4*t*t*w2p))/(2-2*s*s*w2p);
  Ev=(2*t*s*w2d-sqrt(d*d-d*d*s*s*w2d+4*t*t*w2d))/(2-2*s*s*w2d);

  return Ec-Ev-c ;

}

関数f(x,y)の具体的な形

参考までに関数f(x,y)の具体的な形を載せておきます。
イメージ説明
ここで、t=-1,s=0,d=1,eE=0.5として上のグラフは求めましたが、より一般にこれらの値を変化させて等高線の図の変化を考察するのが目的でした。

補足情報

二分法以外のアルゴリズムでの解法がある場合はこちらのご投稿いただけますと幸甚に存じます。

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

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

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

    クリップを取り消します

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

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

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

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

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

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

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

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

    質問の評価を下げる

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

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

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

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

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

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

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

    詳細な説明はこちら

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

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

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

回答 1

checkベストアンサー

0

Evの途中からw2dがw2pになってるのは合ってますか?

二分法の判定のところ・・・f(xmin,y)*f(xmid,y)<0.0 なんですけど、
ここではyは定数として考えることが出来ると思うんです。

つまり、判定は f(x)がxに対して単調であれば、計算による誤差がないという理想的な状況下で
二分法によって正しく求まるとは思いますが、0<=x<=2.5や-2.5<=x<=0の範囲で、
f'(x) は必ずしも正で固定か、負で固定だと言い切れますか?

微分して求めるとこれはわかると思います。確認してもし言い切れないなら、
今回の最適化で二分法を用いること自体が間違っていると思います。

物理学専攻なら、自分で複雑な関数を定義し、最適化の計算方法も細かくプログラムできるようになることはハードルが高すぎると思います。科学技術計算で最近よく使われるPythonを使えるようになるというのはどうでしょうか?信頼性が高く、内部でC言語を呼び出すため、計算もそこそこ速いモジュールがたくさん公開されています。機械学習や大量のデータの整形など、便利なモジュールもあり、これからの科学技術計算のユーザーが使う言語としての将来性はダントツだと思います。

C言語で正しい数値計算を追及しだすと、本業に使うべき時間が奪われることになりそう...
信頼性のある計算方法で出した評価でなければ何の意味もないだろうし...

投稿

  • 回答の評価を上げる

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

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

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

  • 回答の評価を下げる

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

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

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

  • 2016/10/13 09:00

    ご丁寧にありがとうございます。
    Evの式に関して、僕が間違えていました。これを修正したところ、期待していた等高線の図がかけました。

    yが固定されたとき、f'(x)がxによって正になったり負になったりすると、二分法では求まらないということでしょうか。つまり、二分法で求められるのはfが単調増加or単調減少の場合のみということであっていますか?

    自分の所属している研究室ではPtyhonよりもFortranが主流です。Ptyhonでは、上のような計算は比較手容易に可能なモジュールがあるということでしょうか?

    細部まで見ていただきありがとうございます。

    キャンセル

  • 2016/10/13 15:13

    とりあえず解決したみたいでよかったです。

    言葉が足りなかったかもしれません。すいません。
    解を探す範囲として-2.5<=x<=0, 0<=x<=2.5というふうにそれぞれ設定されていますが、二分法で解が少なくともひとつ求まるためには、関数が解を探す範囲[α,β]で連続であって、f(α)*f(β) <=0 あれば十分です。つまり、解を探す範囲として、[α,β]を正しく設定することはとても重要です。ここで、関数が単調であれば、αを限りなく小さい数、βを限りなく大きい数にとることで、解がある場合には正しく求まります。
    つまり、言いたかったことは、-2.5<=x<=0 の範囲で解を探すことをどうやって決定したのか、もしある決定方法をとったのなら、それをほかの可変なパラメータを動かしたときにも対応できるような初期条件の決定方法として一般化することができるのか、ということです。今回の場合、解がひとつ求まればなんでもいいというわけではなさそうです。yが固定でも解として2つ以上のXが対応しうるので...
    となると、正しい初期条件の決定方法を、吟味しながら決めるしかなさそうです。

    私はMATLABについてよく知らなかったのですが、今見てみたら、MATLABはかなり数値計算で余計なことを考えなくでもいいように整備されていると感じました。Pythonを新しく学ぶぐらいなら、すでに使えるMATLABでもいいのかも...
    ただ、速度を気にするならMATLABにしてもやり方を学ぶ必要はありそうです。
    しかし、C言語やFortranとなると途端にハードルが高い感じがします。
    誤差を小さくするために式変形をして、計算順序を考えて、といちいちやるのは厳しすぎますよ

    キャンセル

  • 2016/10/13 15:33

    間違えました。physics303さんが使えるのはMathematicaですね笑
    Mathematicaでもいいんじゃないですか?

    キャンセル

  • 2016/10/13 19:42

    ありがとうございます。よく分かりました。

    MathematicaではFindInstanceという関数をつかってf(x,y)=cをみたす(x,y)を求めることができますが、これはfが複雑になると計算時間が長くなってしまいます。今回、計算するfでは一つ(x,y)をみつけるのに数分かかってしまいます。そのため、CやFortranで数値計算ができないものかなと相談しました。

    もしかしたらMathematicaでもう少し早い方法があるのかもしれませんが、見つけることができませんでした。(私のまわりでMathematicaを使っている人がいなくて・・・)

    キャンセル

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

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

関連した質問

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