/*----------------------------------------*/ /* 熱伝導方程式      */ /*----------------------------------------*/ #include #include #define n 20 /* 一辺あたりのメッシュの数 */ #define T0_amp 10. /* 初期温度分布の中心温度と熱浴の温度の比 */ #define T_edge 1.0e1 /* 境界の熱浴の温度 */ #define Length 1.0e1 /* 矩形の一辺の長さ */ #define kappa 1e0 /* κ */ #define pi 3.141592 /* π */ #define tcal 5. /* 計算時間 */ main(void) { int i,j,k,num; int flag; double factor; double delta_L,dt; double tmp[n+2][2], t; double tout,dtout; FILE *output_file; /* 時間ステップなど */ delta_L = (float)Length/((float)n+1); dt = 0.02*delta_L*delta_L/kappa; factor = dt*kappa/(delta_L*delta_L); output_file=fopen("output1D.dat","w"); /* 出力ファイルオープン*/ /* ---------------- 初期化 及び境界条件 ----------------- */ t=0; tout=0.1; /* 計算を出力する時間 */ dtout=0.1; /* 計算を出力する時間刻み */ /* 境界はメッシュ番号 0 とメッシュ番号 n+1 にとる */ for(i=0; i < n+2 ; i++){ for(k=0; k < 2 ; k++){ tmp[i][k]=T_edge; /* 初期化及び境界の値をT_edgeに */ } } tmp[n/2][0]=T0_amp*T_edge; /* 初期条件設定・真ん中一点が最初高温 */ /* -------------- 初期化 及び境界条件終わり ------------- */ /* -- 初期条件の出力 -- */ for(i=0; i < n+2 ; i++){ fprintf(output_file,"%f %f\n",(float)delta_L*i,tmp[i][0]); } fprintf(output_file,"\n"); fprintf(output_file,"\n"); /* --------------------*/ /* -------- 時間ステップの繰り返しループ--------- */ num=1; while(t<=tcal){ for(i=1; i < n+1 ; i++){ tmp[i][1]=tmp[i][0] + factor*(tmp[i+1][0]+tmp[i-1][0] -2.0*tmp[i][0]); } t=t+dt; num+=1; for(i=1; i < n+1 ; i++){ tmp[i][0]=tmp[i][1]; /* ステップの結果をコピー */ } /* -------- 出力ルーチン -----------*/ if(t>=tout){ /* t >= tout となれば出力 */ printf("number of iteration = %d %f\n",num,t); fprintf(output_file,"#number of iteration = %d %f\n",num,t); for(i=0; i < n+2 ; i++){ fprintf(output_file,"%f %f\n",(float)delta_L*i,(float)tmp[i][0]); } fprintf(output_file,"\n"); fprintf(output_file,"\n"); tout+=dtout; /* toutの更新 */ } /* -------- 出力ルーチン終わり -----------*/ } /* -------- 時間ステップ繰り返しループ 終わり----------- */ }