1次元熱伝導方程式を差分法で解く

質問

モデルは長さL、断面積Aの直方体で
初期温度0度、一端は0度固定でもう一端に熱量Qを与えます。

この熱伝導方程式を差分法(陽解法)で解くプログラム(c言語またはwolfram)を作りたいのですが、上手くいかなくて困っています。

初期条件はT(i,0)=0
境界条件はT(n,j)=0 ←0度固定
T(1,j)=Q*dt/(m*c*dx*A)

差分法の式 T[i][j+1]=T[i][j]+0.5*(T[i+1][j]+T[i-1][j]-2*T[i][j]);

としてfor文で計算させたのですが、いずれ温度が定常状態になるはずが無限に上がり続けてしまいます。

どのようなプログラムを書けば正しく計算できるのでしょうか？

独学で教科書をみながらやっているため気軽に質問できる方が近くにいないためこちらで質問させて頂きました。
なにかアドバイス頂けたら幸いです。
よろしくお願いいたします。

回答No.1

f272 · ベストアンサー

どんなプログラムを書いているの？普通に書けばちゃんと計算できますよ。
#include <stdio.h>
int main(void)
{
int i,j,n=6,maxt=100;
double T[7][101];//T[n+1][maxt+1]
double Q=1,m=1,c=1,A=1,dt=1,dx=1;
for(i=1;i<=n;i++) T[i][0]=0;
printf("tj=%3d",0);
for(i=1;i<=n;i++) printf(" %f",T[i][0]);printf("\n");

for(j=0;j<maxt;j++){
T[n][j+1]=0;
T[1][j+1]=Q*dt/(m*c*dx*A);
for(i=2;i<n;i++) T[i][j+1]=T[i][j]+0.5*(T[i+1][j]+T[i-1][j]-2*T[i][j]);
printf("tj=%3d",j+1);
for(i=1;i<=n;i++) printf(" %f",T[i][j+1]);printf("\n");
}
return 0;
}

補足

回答ありがとうございます。
温度が上がり続けてしまうのは簡単なミスでした。お手数おかけしました。
実行したところ回答者様と同じような結果になりました。
しかしこの結果だと、直方体の一端を熱量qの温度で固定しただけで問題の主旨と合っていないということがわかりました。
また質問内容をあらため質問させていただくと思います。