フィックの第二法則の刻み時間（フォートラン）

質問

第二法則について数値解析を行い、
フォートランによって

dt=1.0
dx=1e-4
d=2e-12
a=d*dt/(dx)**2

do 300 j=0,3600
c(j,0)=0.0
c(j,20)=2.0
do 400 i=1,19
c(j+1,i)=c(j,i)+a*(c(j,i+1)-2.0*c(j,i)+c(j,i-1))
400 continue
300 continue

として一秒ごとに計算し、一時間後までの各時間、各位置の濃度を求めています。
（jは時間、iは位置を表しています。）

このとき、刻み時間t＝１として計算しているのですが、これを0.1秒で計算したいとき、

do 300 j=0,3600 を do 300 j=0,3600,0.1
c(j+1,i) は c(j+0.1,i)

としなくてはいけないのでしょうか？
それとも１のままでよいのでしょうか。
どなたか、どうか教えてください。

ちなみに、上のようにかえてもプログラムが通らないことはわかっています。
聞きたいのは、「刻み時間を変えると濃度計算の中身と計算のステップも変えなくてはいけないのか」ということです。

わかりにくくて申し訳ありません。
どうかお願い致します。

回答No.3

irija_bari · ベストアンサー

ANo.2 です。
例えば、物体の速度 v [m/sec] を例にして考えてみます。
物体の位置を l(t) とすると、
l(t) = l(t-dt) + (v * dt)
（ただし、dt は単位時間（微小時間））
と書くことができますよね。

dt = 1 のときのプログラムの式を、

==== code 1 ====
l(0) = 0
do j = 0, 100
　l(j+1) = l(j) + (v * dt)
continue
==== end of code 1 ====

とすると（インデックスが +1 で時間は +1）、dt = 0.1 のとき（0.1秒刻みで位置を算出するとき）のプログラムの式は、

==== code 2 ====
l(0) = 0
do j = 0, 100, dt
　　l(j/dt+1) = l(j/dt) + (v * dt)
continue
==== end of code 2 ====

となります（インデックスが +1 で時間は +0.1）。
ただ、これでは j/dt が整数にならないとエラーになってしまいますので、以下のように
書き換えるべきではないかと思います。

==== code 3 ====
l(0) = 0
do j = 0, 100/dt
　　l(j+1) = l(j) + (v * dt)
continue
==== end of code 3 ====

"code 2" と "code 3" でやっていることは同じなのですが、イメージとしては、
・"code 2" は dt を元の値の 1/10 にした。
・"code 3" は dt はそのままで速度を 1/10 にした。
という感じです。

私は、"code 3" のように考えるようにしているので、単位が必要と思ったのです。
（あまり必要なかったですが）
う～ん、かえって回りくどい書き方になってしまったかな。。。

お礼

早速回答頂き、どうも有難う御座います。

おぉ！
このプログラムのほうが、私が書いたものよりスッキリしていていいですね！

えぇっと、お書き頂いたプログラムを見る限り、回答としては「時間の刻み幅を変えた場合、濃度計算の時間幅とステップも変えなくてはいけない、ということでよろしいでしょうか・・・？

回答No.2

irija_bari

>a=d*dt/(dx)**2
この式の各変数の単位（例えば、dt なら [sec]）がわかれば、
自ずと決まってくると思います。
もしわからないようでしたら、補足で単位を書き出してみてくれませんか？

お礼

回答いただき、ありがとうございます。

単位ですか・・・。
おそらく、

dは定数で、dxはメートル、dtは秒

です。

どうして単位から決まってくるのでしょうか・・・？
もしよろしければその理由もお教え下さると幸いです。

回答No.1

zettsei

c(*)というのは配列ですよね？
配列の番号に入れるのは整数なので、dtを0.1にしたのならjは1刻みのまま、36000回ループさせるしかありません。

お礼

解答ありがとうございます。

はい、ｄｔを０．１にすると、ご指摘の通り、このプログラムでは上手くいきません。

その問題の解決方法としては、read文でdtを読みとり、if文で、dt=0.1ならばdt=1として36000回ループを行おうと思っておりました。

質問の説明が足りなかったようです。
申し訳ありません・・・。

時間の刻み幅は0.1、1、10の三種類設定しようと思っています。

刻み配列子の番号が整数、という問題は上記の方法でいいとして、ここで問題なのは
「刻み時間を変えると濃度計算の中身と計算のステップも変えなくてはいけないのか」
ということです。

つまり、

read(*,*)dt
read(*,*)n

if(dt.eq.(0.1)) then
k=10
dt=1

else

k=1

end if

do 300 j=0,3600*k,dt
もしくは
do 300 j=0,3600*k

c(j,0)=0.0
c(j,20)=2.0
do 400 i=1,19

c(j+dt,i)=c(j,i)+a*(c(j,i+1)-2.0*c(j,i)+c(j,i-1))
もしくは
c(j+1,i)=c(j,i)+a*(c(j,i+1)-2.0*c(j,i)+c(j,i-1))

400 continue
300 continue

上に二通りのdo文と、二通りの濃度計算がありますが、どちらの計算でやるのか、と言うことです。

回りくどい説明で誠に恐縮なのですが、もう一度お教えいただけたら幸いです。

どうかよろしくお願い致します・・・。

補足

申し訳ありません、｛お礼｝に書いたプログラムは少し間違っていたので・・・。

do 400 i=1,19
ではなく
do 400 i=1,n
の間違いでした。
これでフィックの法則を用いる材料の位置を表しています。

下手な質問でお恥ずかしい限りです・・・。

どうかもう一度お力添えをよろしくお願いします・・・。