MasaYan24Hours -25ページ目

MasaYan24Hours

MasaYanの気まま日記

計算してたら、こんな波打つ数式が出てきた。


MasaYan24Hours-eq

実に滑らかにくねくねしてると思わない???
以下のサイトのやり方を参照して、検索すれば探したいファイルが見つかるかも。

Mac Spotlight検索結果を日本語で条件を指定して絞り込む方法のまとめ / Inforati
http://inforati.jp/apple/mac-tips-techniques/system-hints/how-to-limit-macos-spotlight-search-result-with-parameter-in-japanese.html
ここ2時間くらいずっと計算結果が合わなくて困っていたけどやっと合った。


摂動相互作用を入れた計算結果があるパラメーターの極限で、摂動がないものと一致するようにプログラムを組んだはずなのに、なかなか合わなくて困っていたが、最終的にひらめいた。


2時間かかったけど、少し達成感。


これで、プログラムに不安を抱えずに次に進める。
どこだ。。。どこにバグが。。。

!========== V_G potential ==========
subroutine VGpotential(dVG,il,dr)
include 'head.inc'

real(dp), dimension(mrmax,mrmax), intent(out) :: dVG
integer, intent(in) :: il
real(dp), dimension(mrmax), intent(in) :: dr

real(dp), parameter :: dsigma = 0.1_dp

integer, parameter :: ithe_max = 3000
real(dp), parameter :: dthestep = dpi/dble(ithe_max)

integer :: ir1, ir2
integer :: ithe
real(dp) :: dthe
real(dp) :: dsum

write(55,*) il

write(55,*) dr

dVG = 0._dp

do ir1 = 1, mrmax
do ir2 = 1, mrmax
dsum = 0._dp
do ithe = 1, ithe_max
dthe = dble(ithe) * dthestep
dsum = dsum&
+ cos(dble(il)*dthe) * exp(2d0*dr(ir1)*dr(ir2)*cos(dthe)/dsigma**2)
end do
dVG(ir1,ir2) = 2._dp * dg_G * exp(- (dr(ir1)**2+dr(ir2)**2)/dsigma**2)&
* dsum * dthestep
write(55,*) ir1, ir2, dVG(ir1,ir2)
end do
end do

end subroutine VGpotential


!---
mrmax = 140
il = 0, 1, 2, ...
dr(i) = 0, ..., 16


発見に一時間以上かかった。。。
答え
手計算で整理したときに、exp を二つに分けたんだけど、
実際の計算ではその指数が大きくなりすぎて、オーバーフロー
起こしてたみたい。NaN って形で計算結果がでたけど、
0/0 (ゼロ割るゼロ) 以外にオーバーフローでもその様に
計算結果を返すらしい。
exp の中身が大きくなり過ぎないように、do ループの外の
exp を do ループ内に入れて、足し引きを exp の肩でやれば
解決。
またやってしまいました。


単純なミスによってプログラムが走らなかった。



繰り返すのは愚かなことなので、ここに最初にチェックする基本事項を。



1. 問題があると思われる場所の関係する変数を全て書き出させる (書式を指定せずに)。

  これで、変数の型が実際に正しいものかどうか分かります。

2. 変数は元データ (更新される前のものか) か、コピーデータ (最新のものか) か確認する。


その他は後ほど。。。