C言語による差分法に関して

質問

C言語で非線形偏微分方程式を陽解法による差分法を用いて数値計算しております。

対象となる物体が非常に小さいため、時間刻みを非常に小さくして行っております。
n=40;　 //分割数
d=30e-6;　//物体の長さ
dt=1.0e-15;　//時間刻み幅
loop=1.0e+5;　//ループ回数（最終的にはもっと増やしますが、今はこの程度）

プログラムの詳細な式や値は省略して

double T[?][?] //?のところがわかりません。

for(i=0;i<=40;i++) {
T[i][0]=・・・ //初期条件
}

for(j=0;j<=loop;j++) {
T[0][j]=・・・ //境界条件
T[40][j]=・・・
}

//差分計算
for(j=0;j<=loop;j++) {
for(i=1;i<40;i++) {
T[i][j+1]=・・・・・・・・・・・・・・・・・

}

}

というような流れなのですが、はじめの配列宣言のT[?][?]の？に入れるべき値はなんでしょうか？
分割の数とでループ回数に1をたしたもの（０があるので）かと思いT[41][100001]としたのですが「Segmentation Fault」となってしまいます。一方、T[50][20000]のようなよくわからない数をいれると一応エラーは出ずに動きますが値が変な気がします。（固定したはずの境界条件が変わってしまいます。）
この境界条件が変わってしまうのは配列の宣言に問題があるからのように感じました。

また、もしこの境界条件が変わってしまうのが他の要因にあるとすればどこであると考えられるでしょうか？

この投稿のマルチメディアは削除されているためご覧いただけません。

回答No.4

kmee · ベストアンサー

色々と予想はできますが、プログラム全体を見ないことには、確実なことは言えません。

○配列の大きさ
Cのメモリの使い方はやや難しいので、キーワードを元に参考書や解説サイトを見てください。

double T[?][?] ;
のように、特に何も付けずに宣言したものを「自動変数」と言います。
スタック領域というところに確保されることが多いですが、このスタック領域は、そんなに大きなサイズではありません。
T[41][100001] だと、スタック領域では足りない可能性があります。

これの対策には
・ヒープ領域を使う。具体的には、malloc等を使って確保する
・static(静的変数)を使う。
おそらく、単純なプログラムだろうと思うので、staticで対応できそうな気はします。

○ T[50][20000]
Cでは、配列の範囲外を指定しても、それだけではエラーにはなりません。
その指定した添字から計算されたアドレスをアクセスします。

このアドレスが保護領域で、結果としてエラーになることはあります。
逆に、まったく問題無く動作しているように見えることもあります。

例えば
for(j=0;j<=loop;j++) {
T[0][j]=・・・ //境界条件
T[40][j]=・・・
}
で、j>=20000になる
→配列の範囲外にアクセス
→ そこはループ用の変数j の領域
→ 配列へ代入すると、変数jの領域が書き変わる
→ループ終了条件(j<=loop)をたまたま満してしまった
→ループが終了
と、一見何事もなかったように動作している「ように見える」
ということも起りうることです。

回答No.5

ninoue

回答ではありませんが、少し気になった事があります。
素人が何も分らずに言っているだけかも知れませんが、その場合は申し訳ありません。

>dt=1.0e-15;　//時間刻み幅
偏微分方程式で解くような物理現象は、このように時間刻み幅を小さくしなくても十分ではないでしょうか。
psオーダーでも 1.0e-12ですし、........

各ループ毎の変数セーブ域
前回値と新しい値の変数セーブ域、その他があれば計算出来るのだと思われます。
途中経過も残してアニメ表示等に使われるのでしたら、例えば100ループ、1000ループ等毎にセーブしておいて使うようにされたら良いのではと思われます。

お礼

ご指摘ありがとうございます。
時間刻み幅については私もこんな細かくしたくありません。しかし、上のプログラムに使われているのは偏微分方程式の陽解法という手法なのですが、これは適切にnとdtを選ばないと解が発散してしまうのです。今回はdがマイクロオーダーで小さいため、時間刻みをここまで小さくしないと行けなくなってしまったのです。

セーブに関しては、いいアイディアだと思います。試してみます。

回答No.3

kmee

#1です。訂正です
1e+5=100000 > 20000 ですね。
ということは、#2にある「メモリが割り当てられない」というのがまず疑う点となります

> がT[x][y]のyがループ回数より少ない

ということですが、その際に
n=40;　 //分割数
d=30e-6;　//物体の長さ
dt=1.0e-15;　//時間刻み幅
loop=1.0e+5;　//ループ回数（最終的にはもっと増やしますが、今はこの程度）
というあたりの変数はまったく変えななったのですか?

変えなかったら、同様にSegmentation Faultが 出そうな気がします。よって、プログラム自体が変な可能性があります。
配列に合せて変化調整しているのなら、エラーにはならないかもしれません。

お礼

度々の回答ありがとうございます。

症状をまとめますと
n=40;　 //分割数
d=30e-6;　//物体の長さ
dt=1.0e-15;　//時間刻み幅
loop=1.0e+5;　//ループ回数（最終的にはもっと増やしますが、今はこの程度）
の状態で配列の大きさはT[50][20000]とする。すると、なぜか20000<loopにも関わらず実行ができるが、境界条件が変わってしまう。

上の4つの値は変えずに適切であると考えられる配列の大きさT[41][100001]とするとSegmentation Falt になる。
こんな大きい配列は用意できないということできすかね？最終的には1秒まで、つまり1e+15回ループをしたいのですがそれは無理なんですかね…

回答No.2

f272

Segmentation Faultになってるのは，単に大きな配列を確保しようとしてるからじゃないのかなあ？
後々のことを考えてもmallocとかで動的に領域を確保すべきだと思うよ。

境界条件が変わってしまうのは，これだけじゃわからんが，確かに20000<1.0e+5だからそのせいかもしれない。それ以外の可能性も大いにあると思うけどね。

ついでに変数loopの型は整数型かな？ループ変数を実数型にするのは気持ち悪いし，整数型だとすると1.0e+5を代入するのはやはり気持ち悪い。

回答No.1

kmee

原因と結果が逆ではないか、と予想されます
ループカウンタだけに注目すれば、 T[41][100001] でよいはずですし、
プログラムが正しければ、T[50][20000]としても、余計な部分が使われないだけで、計算結果は同じになるはずだからです。

T[X][Y]のX,Yが常に範囲内か、よく確認してください。
j<=loopのループ内で T[i][j+1] - T[i][j-1] とかやってませんか?

お礼

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

>> j<=loopのループ内で T[i][j+1] - T[i][j-1] とかやってませんか?
そのような間違いはしていませんでした。
ループ回数を1.0e+5としたとき、T[50][20000]でプログラムが実行できますがT[x][y]のyがループ回数より少ないのにエラーがでないのが不思議です。