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

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

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

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

2040閲覧

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

physics303

総合スコア89

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

0クリップ

投稿2016/10/12 11:56

編集2016/10/12 23:55

###前提・実現したいこと
物理学が専門の大学院生です

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

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

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

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

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

C

1include <stdio.h> 2include <math.h> 3double f(double,double); 4 5void main(void){ 6 double x,xmin,xmax,xmid,y,eps; 7 int cnt,cntlimit,i,m; 8 FILE *file; 9 10 eps=0.001; 11 cntlimit=100000; 12 13 file = fopen("data.dat","w"); 14 cnt=0; 15 16 for(i=-500;i<=500;i++){ 17 18 y=2.0+i/100.0; 19 xmax=2.5; 20 xmin=0.0; 21 22 do{ 23 cnt ++; 24 xmid=(xmax+xmin)/2.0; 25 26 if( f(xmin,y)*f(xmid,y)<0.0) xmax=xmid ; 27 else xmin=xmid; 28 29 }while( fabs(xmax-xmin)>eps && cnt < cntlimit ); 30 31 if(cnt==10000) fprintf(file,"It doesnt converge \n"); 32 else if( fabs(f(xmid,y))>0.01 ) fprintf(file,"0\t0\n"); 33 else fprintf(file,"%f\t%f\n",xmid,y); 34 35 36 } 37 38 for(i=-500;i<=500;i++){ 39 40 y=2.0+i/100.0; 41 xmax=-0.0; 42 xmin=-2.5; 43 44 do{ 45 cnt ++; 46 xmid=(xmax+xmin)/2.0; 47 48 if( f(xmin,y)*f(xmid,y)<0.0) xmax=xmid ; 49 else xmin=xmid; 50 51 }while( fabs(xmax-xmin)>eps && cnt < cntlimit ); 52 53 if(cnt==10000) fprintf(file,"It doesnt converge \n"); 54 else if( fabs(f(xmid,y))>0.01 ) fprintf(file,"0\t0\n"); 55 else fprintf(file,"%f\t%f\n",xmid,y); 56 57 58 } 59 60 61 fclose(file); 62} 63 64double f(double kx,double ky){ 65 double c=1.98; 66 double w2p,w2d,eE,t,s,d,Ec,Ev; 67 68 t=-1; 69 s=0; 70 d=1; 71 eE=0.5; 72 73 w2p=3+4*cos(sqrt(3)*(kx+eE)/2.0)*cos(ky/2.0)+2*cos(ky); 74 w2d=3+4*cos(sqrt(3)*(kx-eE)/2.0)*cos(ky/2.0)+2*cos(ky); 75 76 Ec=(2*t*s*w2p+sqrt(d*d-d*d*s*s*w2p+4*t*t*w2p))/(2-2*s*s*w2p); 77 Ev=(2*t*s*w2d-sqrt(d*d-d*d*s*s*w2d+4*t*t*w2d))/(2-2*s*s*w2d); 78 79 return Ec-Ev-c ; 80 81} 82

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

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

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

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

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

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

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

guest

回答1

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/12 16:03

退会済みユーザー

退会済みユーザー

総合スコア0

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

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

physics303

2016/10/13 00:00

ご丁寧にありがとうございます。 Evの式に関して、僕が間違えていました。これを修正したところ、期待していた等高線の図がかけました。 yが固定されたとき、f'(x)がxによって正になったり負になったりすると、二分法では求まらないということでしょうか。つまり、二分法で求められるのはfが単調増加or単調減少の場合のみということであっていますか? 自分の所属している研究室ではPtyhonよりもFortranが主流です。Ptyhonでは、上のような計算は比較手容易に可能なモジュールがあるということでしょうか? 細部まで見ていただきありがとうございます。
退会済みユーザー

退会済みユーザー

2016/10/13 06: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 06:33

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

2016/10/13 10:42

ありがとうございます。よく分かりました。 MathematicaではFindInstanceという関数をつかってf(x,y)=cをみたす(x,y)を求めることができますが、これはfが複雑になると計算時間が長くなってしまいます。今回、計算するfでは一つ(x,y)をみつけるのに数分かかってしまいます。そのため、CやFortranで数値計算ができないものかなと相談しました。 もしかしたらMathematicaでもう少し早い方法があるのかもしれませんが、見つけることができませんでした。(私のまわりでMathematicaを使っている人がいなくて・・・)
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問