「月曜輪講（有田研）」の編集履歴（バックアップ）一覧はこちら
「月曜輪講（有田研）」(2009/07/12 (日) 02:56:57) の最新版変更点
追加された行は緑色になります。
削除された行は赤色になります。
7/11 立川
うまくいきません・・deltakwの初期値を全部1にすると1に就職します・・
&ref(ver.lzh)
7/6 栗田
&ref(attach4_2.zip)
立川君に、zipが開けないと言われたのでもう一度作ってみました。ちなみに、中身のHubbard6,7は定義されてる関数そのものはあまり変わりません。Mとして使っている値が微妙に違うくらいです。
&ref(attach4.zip)
フーリエ変換を改善し、小さな配列で済むようにしました。今までg[N1][N1]だったものがg[n1][n1]で済むようになっています。ちなみにn1=N1/2+1です。それでMを64にしてやってみたけど、やはり+0.25より0.125のほうがいいという結果にｗ
7/5 伊藤
逐次近似でやりました。3回くらいで収束してくれました。初期値の選び方が運が良かったからでしょう。
Green0とGreen(繰りこみ)の比較のためにデータファイルもアップします。
データファイル名は()の中がループをまわした回数です。そのあとが初期値です。
今の状態は、N1が8じゃないとだめです。
あと、今の段階ではグラフは用意していません。
&ref(ver2.zip)
&ref(輪講.pdf)
7/4 伊藤
おそくなってすみません。⑤の繰りこまれたGの計算（ただし、逐次近似）までは作れそうです。
動作確認をし次第アップします。たぶん明日の朝になってしまうかもしれません。
7/1 立川
栗田君ありがとうございました。
Hubbard4-3の件
どういうわけか
if(fp==NULL) return -1
をはずすとうまくうごきます。
Hubbard4-4の件
N1を32
NMAXを1024
NMAXSQRTを32
complex<double> greenw[N1][N1][M], greent[N1][N1][M*2], chiw[N1][N1][M+1], chit[N1][N1][2*M], ab[N1][N1], chi0w[N1][N1][M+1], chi0t[N1][N1][2*M]
の箇所をやめて
chi0w[N1][N1][M+1]だけ確保したらうまくいきました。
やっぱりメモリの使い過ぎみたいですな。
いらなくなったメモリはがんがん解放してかないとだめみたいっす。
ピークの位置の件：
不思議ですね。なんででしょう・・・・。
いまのところomegaをふることとbetaをふることとNをふることと、片っぱしから(double)にキャストすることはやってみたんですが、一向に良くなりません。
なんか突然kx=kyになった瞬間負になるものがあって、それで急激に値が減っているように見えますね。
7/1 伊藤
RPA近似によるギャップ関数の計算の流れについて、
まとめたものを確認のためのせておきます。
&ref(輪講.docx)
①～④までは今まで計算してきた（RPA近似の）χ_0 を使用すれば一瞬でできるはずです。
χ_0はフーリエ変換しないほう（栗田君のchi_0）をとりあえず使おうかと。
chi_s,chi_cもそのまま流用しようと思います。
⑤、⑥のダイソン方程式、ギャップ方程式を解くところは、
ハードルが上がりますが、方程式の形的に逐次代入の方法で求め
るのがよさそうです。
逐次代入がうまくできれば、FLEX近似の方にも応用がきくと思います。
何か意見があったらこちらの方もよろしくお願いします。
↓領域の問題。たぶんプログラムごとに使用できる記憶域が決まっているはずなので、プログラムを分割してしまって、
exeをプログラム中で呼び出すなどすればいけるかも。
6/30 栗田
&ref(class.zip)輪講で使ったファイルを上げておきます。
4_2とあるのが、輪講終了時のファイルです。また、4_3,4_4なども上げたのはいくつか問題となる部分があり、見てもらいたいので上げました。
4_3の問題点：このcppファイルでは現在グラフ描画部分がコメントアウトされていますが、この部分を機能させると、エラーが発生します。謎なのは、この現象がT=0.01などの温度で起こることと、datファイルはちゃんと出力できて、特に発散している様子もなく普通にプロットしてくれそうな値が並んでいることです。ちょっと動かしてみて、エラーが出るかとか確かめてみてください。報告お待ちしてます。
4_4の問題点：(4.14)に基づいてchi_0を求めてくれるのですが、微妙にピークが出るのはkx=kyのときではなく、|kx-ky|=dkのときです。グラフ描画のコメントアウトを外してプロットしてもわかりますし、datファイルを見てもわかると思います。chi_0の求め方は問題ないはずなんですが、1週間あるのでどなたかチェックしてみてください。
一般的な問題点：N1=32とか、M=64とかやるとエラーが出ます。多分complex<double> greenw[N1][N1][M], greent[N1][N1][M*2], chiw[N1][N1][M+1], chit[N1][N1][2*M], chi0w[N1][N1][M+1], chi0t[N1][N1][2*M];の領域を確保するだけで、コンピューターの限界なんじゃないかと思います。T=0.01でやることにすると、M=64でも足りないくらいなのにこまったものです…。
僕のブラウザだと、月曜輪講というタイトルの上にものすごい空白があったので、立川君の6/22のコメントを別ページに移動しました。wikiモードのときにプログラムコードをそのまま書くと、C++のコードをwikiが変なふうに解釈するので、Cのコードとかをそのまま書くときは、シンプルテキストモードの新規ページを作ってそこにリンクしてください。
6/29 立川
グラフ機能が付きました
&ref(attach2-graph(2).lzh)
6/29 栗田
&ref(attach3.zip) パデ近似のアルゴリズムをあげて起きます。遅くなってすいません。
6/28 栗田
&ref(attach2.zip) フーリエ変換をしてχ0を求めるプログラムを上げておきます。実際のχ0(4.14)式と一緒に表示するようになっています。日曜に目を通してもらえると助かります。日曜中にコメントがあればなお嬉しいです。あと、前回の輪講が終わった時点での定数（化学ポテンシャルとか温度とか格子定数とか）を教えてもらえると嬉しいのですが・・・。
6/22 立川
[[gnuplotをC++で操作する方法>メモ1]]の詳細を説明いたします。
6/20 栗田
C++のプログラムを書いたので上げときます。&ref(attach1.zip) &ref(attach1_Readme.txt)
立川君がメールでくれたプログラムで、少し気になったところがあったので、自分だったらこうするという発想に基づいて作ったものです。気になったところは以下のところです。多分僕が何か勘違いをしている部分もあると思いますが。
グリーン関数の引数wについて。複素変数にしました。実数のwで複雑な連立方程式（エリアシュベルグみたいなやつ）を解こうとすると、パラメータによってG0が発散しそうでちょっと怖いです。そもそも、エリアシュベルグみたいな方程式が実周波数空間で成り立つかはよくわからないので…。具体的には、complex.hというヘッダファイルをincludeすれば使えます。
chi_0の導出について。立川君が作ったやつの導出がわからなかったんだけど、あれは畳み込みの性質を使ってるんですか？(4.14)の式が微妙に畳み込みと違っているのが気になるんですが…。もっと気になるのは、畳み込みだとしてchi0=G0*G0が成り立つのは(r,τ)空間だと思います。フーリエ変換周辺はプログラムが難しくて断定はできないけど、松原周波数をτになおしてないように見えます（つまりchi0(r,iw)=G0(r,iw)*G0(r,iw)をやっているように見える）。いずれにせよ(4.14)の下の式（フェルミ関数とかが出てくるやつ）を使ったほうが早いと思います。
逆フーリエ変換について。伊藤君が上げてくれたフーリエ変換のプログラムだけど、example1.cではcdft2d(-1)（つまり逆フーリエ変換）を呼んだ後で、要素数（2次元ならN*N)で割っていますが、これが抜けていると思います。
物理定数の定義について。格子定数があったほうがわかりやすい気がします。Hubbard3.cppの16行目のcos(kx)という表現から、格子定数は1を採用しているように見えますが、こうすると少し不都合が起こると思います。というのは、Nを有限にして周期的境界条件を考える場合kは離散的なベクトルになりますが、このときの刻み幅は格子定数およびサイト数に依存します（アシュクロフト22・26など）。つまり、N,a,dkは独立には選ばないほうが無難だと思います。kについて連続変数みたいなものなので、刻み幅は大して問題になりませんが、考えるべきkの区間の広さもも格子定数によると思うので注意したほうがいいと思います。
6/20 栗田
立川君が上げたと思われるsciファイルが埋もれてるので上げときます。
&ref(test.sci)
6/19 立川
まえ要望があったので、Ｃ＋＋でgnuplotをコントロールするコードを書いておきます。
ソースとおなじフォルダのpgnuplot.exeというファイルをおいてください
wgnuplot.exeとは違うものです。
ちなみに私がお送りしたコードでは、
datファイルに計算結果をかきこむ
pltファイルを書かせる(gnuplotで読み込むとグラフをそのとおりに書いてくれる）
pltファイルをpgnuplotで実行というスタイルをとっています。
pgnuplot.exeをC++でコントロールする方法のテスト
FILE *fp=_popen("pgnuplot.exe","w");
if(fp==NULL)
return -1;
fputs("plot sin(x)\n",fp);
fflush(fp);
cin.get();
_pclose(fp);
6/16 立川
こんばんは。岩澤君が遊びに来ていますｗｗｗｗｗ
２次元フーリエ変換まで実行して、χ_sの波数依存性まで計算いたしました。
scilabとC++(Ｃを用いた範囲で解決しています）
ソースコードが見やすくないと思うので読めないかもしれません・・・
そのうち調整します。
6/13伊藤
2次元FFT（http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html にあるもののうち、二次元 複素離散 Fourier 変換、のみを使う場合）
の使用方法について。
&ref(fft2dc.h)
を使います。ヘッダーファイルとして読み込んでおいてください。
2次元配列（N＊Nを想定）のFFTのやりかた↓
&ref(example1.c)
ダウンロードしてから見てください。
これが、使いにくかったら、栗田くんが前に使ったやつ
http://www5.airnet.ne.jp/tomy/cpro/sslib8.htm
のやつの方が楽かも（特に配列を用意する必要はなかったはず。）
6/09 栗田
五月祭の磁性体班の情報交換用に作ったwikiに有田研の輪講ページを作りましたｗ
次回からは数値計算っぽいこともやるようなので、担当になった人は事前にプログラムとかアップしてくれると助かります。ところで、C言語を使うことになりそうですが、ぶっちゃけscilabでもできるような・・・。まあCでもいいですが。Cでやるならグラフ化とかの説明よろしくお願いします。
7/12 栗田
少しプログラムを見てみました。個人的にDeltaはもう少し小さい数、complex<double> (0,0.000000001)くらいのほうがいいと思います。あと、あがっているプログラムでは方程式を固有値問題と考えている感じですが、実際に固有値がどのくらいになるか計算した方がいいと思います。つまりある程度計算が進んだら、standardで割らずにdeltaknewがdeltakwの何倍か見て、その固有値に応じてTを変えていくのが、先生が言ってたことかと。ただ、松原周波数が少ないので、望むような結果がでないときアルゴリズムが悪いのか周波数が少ないのかよくわからないよね・・・。あまり役に立たないコメントですいません。あと、
http://www.hiroshima-cu.ac.jp/japanese/IPC/hunet99/sun/WorkShop/ja/html_docs/c-plusplus/c++_lrm/Chapter3.doc.html
を見ると、conjとabsは複素数に対して定義されているので、良かったら使うといいかも。
7/11 立川
うまくいきません・・deltakwの初期値を全部1にすると1に就職します・・
&ref(ver.lzh)
7/6 栗田
&ref(attach4_2.zip)
立川君に、zipが開けないと言われたのでもう一度作ってみました。ちなみに、中身のHubbard6,7は定義されてる関数そのものはあまり変わりません。Mとして使っている値が微妙に違うくらいです。
&ref(attach4.zip)
フーリエ変換を改善し、小さな配列で済むようにしました。今までg[N1][N1]だったものがg[n1][n1]で済むようになっています。ちなみにn1=N1/2+1です。それでMを64にしてやってみたけど、やはり+0.25より0.125のほうがいいという結果にｗ
7/5 伊藤
逐次近似でやりました。3回くらいで収束してくれました。初期値の選び方が運が良かったからでしょう。
Green0とGreen(繰りこみ)の比較のためにデータファイルもアップします。
データファイル名は()の中がループをまわした回数です。そのあとが初期値です。
今の状態は、N1が8じゃないとだめです。
あと、今の段階ではグラフは用意していません。
&ref(ver2.zip)
&ref(輪講.pdf)
7/4 伊藤
おそくなってすみません。⑤の繰りこまれたGの計算（ただし、逐次近似）までは作れそうです。
動作確認をし次第アップします。たぶん明日の朝になってしまうかもしれません。
7/1 立川
栗田君ありがとうございました。
Hubbard4-3の件
どういうわけか
if(fp==NULL) return -1
をはずすとうまくうごきます。
Hubbard4-4の件
N1を32
NMAXを1024
NMAXSQRTを32
complex<double> greenw[N1][N1][M], greent[N1][N1][M*2], chiw[N1][N1][M+1], chit[N1][N1][2*M], ab[N1][N1], chi0w[N1][N1][M+1], chi0t[N1][N1][2*M]
の箇所をやめて
chi0w[N1][N1][M+1]だけ確保したらうまくいきました。
やっぱりメモリの使い過ぎみたいですな。
いらなくなったメモリはがんがん解放してかないとだめみたいっす。
ピークの位置の件：
不思議ですね。なんででしょう・・・・。
いまのところomegaをふることとbetaをふることとNをふることと、片っぱしから(double)にキャストすることはやってみたんですが、一向に良くなりません。
なんか突然kx=kyになった瞬間負になるものがあって、それで急激に値が減っているように見えますね。
7/1 伊藤
RPA近似によるギャップ関数の計算の流れについて、
まとめたものを確認のためのせておきます。
&ref(輪講.docx)
①～④までは今まで計算してきた（RPA近似の）χ_0 を使用すれば一瞬でできるはずです。
χ_0はフーリエ変換しないほう（栗田君のchi_0）をとりあえず使おうかと。
chi_s,chi_cもそのまま流用しようと思います。
⑤、⑥のダイソン方程式、ギャップ方程式を解くところは、
ハードルが上がりますが、方程式の形的に逐次代入の方法で求め
るのがよさそうです。
逐次代入がうまくできれば、FLEX近似の方にも応用がきくと思います。
何か意見があったらこちらの方もよろしくお願いします。
↓領域の問題。たぶんプログラムごとに使用できる記憶域が決まっているはずなので、プログラムを分割してしまって、
exeをプログラム中で呼び出すなどすればいけるかも。
6/30 栗田
&ref(class.zip)輪講で使ったファイルを上げておきます。
4_2とあるのが、輪講終了時のファイルです。また、4_3,4_4なども上げたのはいくつか問題となる部分があり、見てもらいたいので上げました。
4_3の問題点：このcppファイルでは現在グラフ描画部分がコメントアウトされていますが、この部分を機能させると、エラーが発生します。謎なのは、この現象がT=0.01などの温度で起こることと、datファイルはちゃんと出力できて、特に発散している様子もなく普通にプロットしてくれそうな値が並んでいることです。ちょっと動かしてみて、エラーが出るかとか確かめてみてください。報告お待ちしてます。
4_4の問題点：(4.14)に基づいてchi_0を求めてくれるのですが、微妙にピークが出るのはkx=kyのときではなく、|kx-ky|=dkのときです。グラフ描画のコメントアウトを外してプロットしてもわかりますし、datファイルを見てもわかると思います。chi_0の求め方は問題ないはずなんですが、1週間あるのでどなたかチェックしてみてください。
一般的な問題点：N1=32とか、M=64とかやるとエラーが出ます。多分complex<double> greenw[N1][N1][M], greent[N1][N1][M*2], chiw[N1][N1][M+1], chit[N1][N1][2*M], chi0w[N1][N1][M+1], chi0t[N1][N1][2*M];の領域を確保するだけで、コンピューターの限界なんじゃないかと思います。T=0.01でやることにすると、M=64でも足りないくらいなのにこまったものです…。
僕のブラウザだと、月曜輪講というタイトルの上にものすごい空白があったので、立川君の6/22のコメントを別ページに移動しました。wikiモードのときにプログラムコードをそのまま書くと、C++のコードをwikiが変なふうに解釈するので、Cのコードとかをそのまま書くときは、シンプルテキストモードの新規ページを作ってそこにリンクしてください。
6/29 立川
グラフ機能が付きました
&ref(attach2-graph(2).lzh)
6/29 栗田
&ref(attach3.zip) パデ近似のアルゴリズムをあげて起きます。遅くなってすいません。
6/28 栗田
&ref(attach2.zip) フーリエ変換をしてχ0を求めるプログラムを上げておきます。実際のχ0(4.14)式と一緒に表示するようになっています。日曜に目を通してもらえると助かります。日曜中にコメントがあればなお嬉しいです。あと、前回の輪講が終わった時点での定数（化学ポテンシャルとか温度とか格子定数とか）を教えてもらえると嬉しいのですが・・・。
6/22 立川
[[gnuplotをC++で操作する方法>メモ1]]の詳細を説明いたします。
6/20 栗田
C++のプログラムを書いたので上げときます。&ref(attach1.zip) &ref(attach1_Readme.txt)
立川君がメールでくれたプログラムで、少し気になったところがあったので、自分だったらこうするという発想に基づいて作ったものです。気になったところは以下のところです。多分僕が何か勘違いをしている部分もあると思いますが。
グリーン関数の引数wについて。複素変数にしました。実数のwで複雑な連立方程式（エリアシュベルグみたいなやつ）を解こうとすると、パラメータによってG0が発散しそうでちょっと怖いです。そもそも、エリアシュベルグみたいな方程式が実周波数空間で成り立つかはよくわからないので…。具体的には、complex.hというヘッダファイルをincludeすれば使えます。
chi_0の導出について。立川君が作ったやつの導出がわからなかったんだけど、あれは畳み込みの性質を使ってるんですか？(4.14)の式が微妙に畳み込みと違っているのが気になるんですが…。もっと気になるのは、畳み込みだとしてchi0=G0*G0が成り立つのは(r,τ)空間だと思います。フーリエ変換周辺はプログラムが難しくて断定はできないけど、松原周波数をτになおしてないように見えます（つまりchi0(r,iw)=G0(r,iw)*G0(r,iw)をやっているように見える）。いずれにせよ(4.14)の下の式（フェルミ関数とかが出てくるやつ）を使ったほうが早いと思います。
逆フーリエ変換について。伊藤君が上げてくれたフーリエ変換のプログラムだけど、example1.cではcdft2d(-1)（つまり逆フーリエ変換）を呼んだ後で、要素数（2次元ならN*N)で割っていますが、これが抜けていると思います。
物理定数の定義について。格子定数があったほうがわかりやすい気がします。Hubbard3.cppの16行目のcos(kx)という表現から、格子定数は1を採用しているように見えますが、こうすると少し不都合が起こると思います。というのは、Nを有限にして周期的境界条件を考える場合kは離散的なベクトルになりますが、このときの刻み幅は格子定数およびサイト数に依存します（アシュクロフト22・26など）。つまり、N,a,dkは独立には選ばないほうが無難だと思います。kについて連続変数みたいなものなので、刻み幅は大して問題になりませんが、考えるべきkの区間の広さもも格子定数によると思うので注意したほうがいいと思います。
6/20 栗田
立川君が上げたと思われるsciファイルが埋もれてるので上げときます。
&ref(test.sci)
6/19 立川
まえ要望があったので、Ｃ＋＋でgnuplotをコントロールするコードを書いておきます。
ソースとおなじフォルダのpgnuplot.exeというファイルをおいてください
wgnuplot.exeとは違うものです。
ちなみに私がお送りしたコードでは、
datファイルに計算結果をかきこむ
pltファイルを書かせる(gnuplotで読み込むとグラフをそのとおりに書いてくれる）
pltファイルをpgnuplotで実行というスタイルをとっています。
pgnuplot.exeをC++でコントロールする方法のテスト
FILE *fp=_popen("pgnuplot.exe","w");
if(fp==NULL)
return -1;
fputs("plot sin(x)\n",fp);
fflush(fp);
cin.get();
_pclose(fp);
6/16 立川
こんばんは。岩澤君が遊びに来ていますｗｗｗｗｗ
２次元フーリエ変換まで実行して、χ_sの波数依存性まで計算いたしました。
scilabとC++(Ｃを用いた範囲で解決しています）
ソースコードが見やすくないと思うので読めないかもしれません・・・
そのうち調整します。
6/13伊藤
2次元FFT（http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html にあるもののうち、二次元 複素離散 Fourier 変換、のみを使う場合）
の使用方法について。
&ref(fft2dc.h)
を使います。ヘッダーファイルとして読み込んでおいてください。
2次元配列（N＊Nを想定）のFFTのやりかた↓
&ref(example1.c)
ダウンロードしてから見てください。
これが、使いにくかったら、栗田くんが前に使ったやつ
http://www5.airnet.ne.jp/tomy/cpro/sslib8.htm
のやつの方が楽かも（特に配列を用意する必要はなかったはず。）
6/09 栗田
五月祭の磁性体班の情報交換用に作ったwikiに有田研の輪講ページを作りましたｗ
次回からは数値計算っぽいこともやるようなので、担当になった人は事前にプログラムとかアップしてくれると助かります。ところで、C言語を使うことになりそうですが、ぶっちゃけscilabでもできるような・・・。まあCでもいいですが。Cでやるならグラフ化とかの説明よろしくお願いします。
