/*----------------------------------------*/ /* 太鼓の膜の振動     */ /*----------------------------------------*/ #include #include #define n 20 /* 一辺あたりのメッシュの数 */ #define max_it 1000 /* 最大時間ステップ数    */ #define u0_amp 1.0e0 /* 初期条件の振幅 */ #define Length 1.0e1 /* 矩形の一辺の長さ */ #define cs 1e0 /* 波の伝播速度 */ #define pi 3.141592 /* π */ main(void) { int i,j,k,num; int flag; double factor; double delta_L,dt; double u[n+2][n+2][3], t; /* u の最後の添え字0,1,2はテキストのk-1,k,k+1に対応する */ FILE *output_file; /* 時間ステップなど */ delta_L = (float)Length/((float)n+1); dt = 0.1*delta_L/cs; factor=pow(cs*dt/delta_L,2); printf("hoge1\n"); output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ /* ---------------- 初期化 及び境界条件 ----------------- */ t=0; /* 境界はメッシュ番号 0 とメッシュ番号 n+1 にとる */ for(i=0; i < n+2 ; i++){ for(j=0; j < n+2 ; j++){ for(k=0; k < 3 ; k++){ u[i][j][k]=0.0; /* 初期化及び境界の値を0に */ } } } for(i=1; i < n+1 ; i++){ for(j=1; j < n+1 ; j++){ u[i][j][1]=u0_amp*sin(pi*(float)i/((float)n+1)) *sin(pi*(float)j/((float)n+1)); /* 初期条件設定 */ } } /* この段階でt=0はk=1に対応することに注意 */ /* -------------- 初期化 及び境界条件終わり ------------- */ /* -- 初期条件の出力 -- */ for(i=0; i < n+2 ; i++){ for(j=0; j < n+2 ; j++){ fprintf(output_file,"%f %f %f\n",(float)delta_L*i,(float)delta_L*j,u[i][j][1]); } fprintf(output_file,"\n"); } /* --------------------*/ /* --------1ステップ目だけは微分の境界条件を使うので別に -----*/ for(i=1; i < n+1 ; i++){ for(j=1; j < n+1 ; j++){ u[i][j][2]=u[i][j][1] + 5e-1*factor* (u[i+1][j][1]+u[i-1][j][1]+u[i][j+1][1]+u[i][j-1][1] -4.0*u[i][j][1]); } } t=t+dt; /* -------1ステップ目終わり----------------------------------*/ /* -------- 時間ステップの繰り返しループ(2ステップ目以降) --------- */ for(num=2; num < max_it; num++){ for(i=1; i < n+1 ; i++){ for(j=1; j < n+1 ; j++){ u[i][j][0]=u[i][j][1]; u[i][j][1]=u[i][j][2]; /* 前のステップの結果をコピー */ } } for(i=1; i < n+1 ; i++){ for(j=1; j < n+1 ; j++){ u[i][j][2]=2.0*u[i][j][1]-u[i][j][0] + factor* (u[i+1][j][1]+u[i-1][j][1]+u[i][j+1][1]+u[i][j-1][1] -4.0*u[i][j][1]); } } t=t+dt; /* -------- 出力ルーチン -----------*/ if(num%10==0){ /* 10 step ごとに出力 */ printf("number of iteration = %d\n",num); for(i=0; i < n+2 ; i++){ for(j=0; j < n+2 ; j++){ fprintf(output_file,"%f %f %f\n",(float)delta_L*i,(float)delta_L*j,(float)u[i][j][1]); } fprintf(output_file,"\n"); } } /* -------- 出力ルーチン終わり -----------*/ } /* -------- 時間ステップ繰り返しループ 終わり----------- */ }