Q&A
提示のコードはどういう動作をするんでしょうか
前提・実現したいこと
砂紋のパターン形成のシミュレーションを行いたい。
発生している問題・エラーメッセージ
上記のシミュレーションを行うプログラムを作成しています。
プログラム上で周期境界を再現したいのですが、やり方がどうしても分かりません。
簡単に言うと、x=50,y=50の幅の平面で2次元配列h[50][50]を扱う際に
[-1]など中身が負の数になってしまうとき、[-1]→[49]として[-1]をもう一方の端として扱うようにしたいです。
C言語に慣れておらず拙い説明で申し訳ありませんが、ご回答いただけると嬉しいです。
該当のソースコード
#include<stdio.h> #include<stdlib.h> #include<time.h> int main(void) { int x_scale, y_scale, step,ir,jr; x_scale = 50; //x方向長さ y_scale = 50; //y方向長さ step = 1; //計算回数 double yuragi, q, b, D, side, slant,v,w; yuragi = 0.3; // q = 0.1; //跳躍により移動する砂の量 b = 1.0; //跳躍距離のパラメータ D = 0.1; //転がりにより移動する砂の量のパラメータ(転がり出る砂の量=Dh) side = 2.0; //転がり計算の時に使うパラメータ slant = 1.0; //同上 srand(time(NULL)); int i,j,s, t, L, k, L_0;//L=L_0 + bh (跳躍距離) L_0 = 3; //風速に対応するパラメータ double h[50][50]; //[x,y]における砂の高さ(積載量) FILE*fp; fp=fopen("output.dat", "w"); for (j = 0;j < y_scale;j++) { for (i = 0;i < x_scale;i++) { h[i][j] = 2.85 + (yuragi * rand() / (double)RAND_MAX); //各座標のhの初期値を2.85~3.15の間で乱数で取る } } for (t = 1; t <= step; t++) { for (k = 1; k < 2500; k++) //跳躍計算 { v = (int)(50 * rand() / (double)RAND_MAX); w = (int)(50 * rand() / (double)RAND_MAX); ir = v; jr = w; // if (h[ir][jr] > 0) { L = L_0 + (b * (int)h[(ir)][(jr)]); h[(ir)][(jr)] -= q; h[(ir + L)][(jr)] += q; } else { continue; } } for (j = 1;j < y_scale - 1;j++) { for (i = 1;i < x_scale - 1;i++) { if (h[i][j] > h[(i - 1)][(j)] && h[i][j] > h[(i + 1)][(j)] && h[i][j] > h[(i)][(j - 1)] && h[i][j] > h[i][(j + 1)] && h[i][j] > h[(i - 1)][(j - 1)] && h[i][j] > h[i - 1][(j + 1)] && h[i][j] > h[(i + 1)][(j - 1)] && h[i][j] > h[i + 1][(j + 1)]) { h[(i - 1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i + 1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i - 1)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i - 1)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i + 1)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i + 1)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[i][j] -= D * h[i][j]; } else { continue; } } } for (j = 0;j < y_scale;j++) { for(i = 0; i < x_scale; i++ ){ printf("%f\n",h[i][j]); fprintf(fp, "%f\n",h[i][j]); } } } fclose(fp); }
変更後のソースコード
#include<stdio.h> #include<stdlib.h> #include<time.h> int main(void) { int x_scale, y_scale, step; x_scale = 50; //x方向長さ y_scale = 50; //y方向長さ step = 1; //計算回数 double yuragi, q, b, D, side, slant,v,w; yuragi = 0.3; // q = 0.1; //跳躍により移動する砂の量 b = 1.0; //跳躍距離のパラメータ D = 0.1; //転がりにより移動する砂の量のパラメータ(転がり出る砂の量=Dh) side = 2.0; //転がり計算の時に使うパラメータ slant = 1.0; //同上 srand(time(NULL)); int x,y,i,j,i_1,i_2,j_1,j_2,i_0,s, t, L, k, L_0;//L=L_0 + bh (跳躍距離) L_0 = 3; //風速に対応するパラメータ i = (x_scale + (x % x_scale))%(x_scale); j = (y_scale + (y % y_scale))%(y_scale); i_1 = (x_scale + ((x - 1) % x_scale))%(x_scale); i_2 = (x_scale + ((x + 1) % x_scale))%(x_scale); j_1 = (y_scale + ((y - 1) % y_scale))%(y_scale); j_2 = (y_scale + ((y + 1) % y_scale))%(y_scale); i_0 = (x_scale + ((x + L) % x_scale))%(x_scale); double h[50][50]; //[x,y]における砂の高さ(積載量) FILE*fp; fp=fopen("output.dat", "w"); for (y = 0;y < y_scale;y++) { for (x = 0;x < x_scale;x++) { h[x][y] = 2.85 + (yuragi * rand() / (double)RAND_MAX); //各座標のhの初期値を2.85~3.15の間で乱数で取る } } for (t = 1; t <= step; t++) { for (k = 1; k < 2500; k++) //跳躍計算 { v = (int)(50 * rand() / (double)RAND_MAX); w = (int)(50 * rand() / (double)RAND_MAX); x = v; y = w; L = L_0 + (b * (int)h[(i)][(j)]); if (h[i][j] > 0) { h[(i)][(j)] -= q; h[(i_0)][(j)] += q; } else { continue; } } for (y = 0;y < y_scale ;y++) { for (x = 0;x < x_scale ;x++) { if (h[i][j] > h[(i_1)][(j)] && h[i][j] > h[(i_2)][(j)] && h[i][j] > h[(i)][(j_1)] && h[i][j] > h[i][(j_2)] && h[i][j] > h[(i_1)][(j_1)] && h[i][j] > h[i_1][(j_2)] && h[i][j] > h[(i_2)][(j_1)] && h[i][j] > h[i_2][(j_2)]) { h[(i_1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i_2)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant)); h[(i_1)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i_1)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i_2)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[(i_2)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant)); h[i][j] -= D * h[i][j]; } else { continue; } } } for (y = 0;y < y_scale;y++) { for(x = 0; x < x_scale; x++ ){ printf("%f\n",h[i][j]); fprintf(fp, "%f\n",h[i][j]); } } } fclose(fp); }
出力
-0.039956
-0.039956
-0.039956
-0.039956
↑これがずっと続きます
補足情報（FW/ツールのバージョンなど）
Visual Studio Codeで作成しています。
動作はするのですが、出力されたファイルに何も書かれていない、といった感じです。
あと、コードは"<code>"または"<コード>"のボタンがあるはずなので、それを押して出てくるヤツの中に書きましょう。今のままでは読みづらいです。
質問は編集できるので編集しましょう。
> 動作はするのですが、出力されたファイルに何も書かれていない
計算うんぬんでは無く、ファイルが空になってしまうという問題なのでしょうか。
>>BeatStarさん
すいません、編集しました。
majuponiさんの回答を参考に、プログラムを編集しました。
数値は出力されるようになったのですが、出力された数値がすべて同じで、しかも負の値になってしまいます。(0以上を想定していた）
おそらく乱数の使い方に誤りがあるとおもうのですが、よくわかりません。。
rand() による乱数の求め方を間違っていないでしょうか。
rand() のドキュメントや、 rand() を含む計算の結果を printf で確認しては如何でしょう。
rand()を含む計算の結果を確認したところ、やはりすべて同じ数値になりました。
> すべて同じ数値
・・・もう一声・・・といいますか、つまり乱数を利用する計算が間違っている、rand() の利用方法が間違っているとは言えないでしょうか?
[Cライブラリ関数-rand（）](https://www.tutorialspoint.com/c_standard_library/c_function_rand.htm )
回答2件
0
周期境界とやらは存じませんが、今回の場合、そのアルゴリズムについての質問では無い為回答します。
...[-1]など中身が負の数になってしまうとき、[-1]→[49]として[-1]...
つまり、『負の数でアクセスする場合、強制的に最大値にする』的な意味ですよね?
それなら単純にそう実装すればいいだけでは。
そもそもプログラミングというのは『こう書けばいい』というものではありません。
言うなれば『現実世界のシミュレーション』です。
たとえば、この周期境界のやつを実際に手作業でやってみてください。
そして、今回の質問である、『負の数なら~』も同様にやってみてください。
そうすると、『要素数に相当する変数が負の数の場合』-> 『要素数の変数を最大値にする』となるはずです。
たとえば、
int pos = -1; int desk[10][10]; int max = 10; // 各要素の最大値 // posが-1なら if( pos == -1 ){ // 強制的に最大値にする pos = max; } // 取得する int n = desk[0][pos];
的な感じで。
ただし、負の数っていうのは -1 だけじゃなくて -2, -3, -4, ... とあります。
それをすべてやるのは非効率。何か規則性があるはずです。
数直線を思い浮かべてください。実際にノートに書いてもいいです。
とにかく、負の数を出来る限りチェックつけましょう。-1, -2, -3, -4, ... と。
で、そうすると『常に0より小さい』ことがわかりますね。そう、『0より小さいかどうか』で判断できますね。0は…含みませんね。
ということは、 pos＜0 の条件式を満たす場合が負の数ですね。
さあ、これでわかりましたね。後は組み込んで、デバッグするだけです。
投稿2021/11/06 11:52
総合スコア4962
あなたの回答
tips
プレビュー
