c言語のプログラムを実行すると強制終了してしまう

Question

前提・実現したいこと

win7のgcc環境において、以下のプログラムでルンゲクッタ法による数値計算を行おうとしています。gccでのデバッグは問題ないのですが、プログラムを実行すると強制終了を起こしてしまいます。原因が特定できていませんので、修正方法をご教授頂きたいです。

#include <stdio.h>
#include <math.h>

double f1(double t,double u,double v);
double f2(double t,double u,double v);

int main(void) {
    FILE *output;
    output = fopen("output.dat","w");
    
    double t = 0.0;
    double dt = 0.01;
    double delta = 0.1;
    double r[3000][3000];
    double u = 1.0, w = 1.0;//初期値
    double a=0.03, b=0.10;//論文での値
    double k0[2][2], k1[2][2], k2[2][2], k3[2][2];//計算した値を入れる配列

    //離散化を構成する
    int i;
    int n,s,p;//mpとmqの添え字
    double mp,mq;
    double u2[3000][3000];
    double w2[3000][3000];
    //double u2[0][0]=0.0;//初期化
    //double w2[0][0]=0.0;//初期化
    u2[s][n]=1.0;
    w2[s][n]=1.0;

    
    for(n=0;n<=3000;n++){
        for(s=0;s<=3000;s++){
            for(n=0;n<=3000;n++){
                u2[s+1][n] = (mp*u2[s][n] + delta*(2*mp*u2[s][n]*w2[s+1][n]+a) )/1 + delta*(w2[s+1][n]*w2[s+1][n]+a+1);
                w2[s+1][n] = (mq*w2[s][n] + delta*(mp*u2[s][n]*(mq*w2[s][n]*w2[s][n]+1)+b) )/1 + delta*(w2[s+1][n]*w2[s+1][n]+a+1);
            }
        }
        mp = 0.5*(u2[s][n+1]+ u2[s][n-1]);
        mq = 0.5*(u2[s][n+1]+ u2[s][n-1]);

    for(t = 0.0; t <=3000; t+=dt) {//ルンゲクッタ法で計算
        k0[0][0]=dt*f1(t,u,w);
        k0[1][1]=dt*f2(t,u,w);
        
        k1[0][0]=dt*f1(t+dt/2.0,u2[s+1][n]+k0[0][0]/2.0,w2[s+1][n]+k0[1][1]/2.0);
        k1[1][1]=dt*f1(t+dt/2.0,u2[s+1][n]+k0[0][0]/2.0,w2[s+1][n]+k0[1][1]/2.0);
        
        k2[0][0]=dt*f1(t+dt/2.0,u2[s+1][n]+k1[0][0]/2.0,w2[s+1][n]+k1[1][1]/2.0);
        k2[1][1]=dt*f1(t+dt/2.0,u2[s+1][n]+k1[0][0]/2.0,w2[s+1][n]+k1[1][1]/2.0);
        
        k3[0][0]=dt*f1(t+dt,u2[s+1][n]+k2[0][0],w2[s+1][n]+k2[1][1]);
        k3[1][1]=dt*f1(t+dt,u2[s+1][n]+k2[0][0],w2[s+1][n]+k2[1][1]);
        
        u2[s+1][n] = u2[s][n] + (k0[0][0]+ k1[0][0]*2 + k2[0][0]*2 + k3[0][0])/6.0;
        w2[s+1][n] = w2[s][n] + (k0[1][1]+ k1[1][1]*2 + k2[1][1]*2 + k3[1][1])/6.0;

        fprintf(output,"%lf, %lf, %lf\n", u2[s+1][n] , w2[s+1][n], t);
        printf("%lf , %lf , %lf\n",u2[s+1][n] , w2[s+1][n], t);
        }
      }
    fclose(output);
    return 0;
}


double f1(double t,double u,double w){
    int n,s;
    double r[3000][3000];
    double a=0.04;
    double u2[3000][3000];
    double w2[3000][3000];
    r[s+1][n] = -u2[s+1][n]*w2[s+1][n]*w2[s+1][n]+a*(1-u2[s+1][n]);
    return 0;
}

double f2(double t,double u,double w){
    int n,s;
    double r[3000][3000];
    double b=0.11;
    double u2[3000][3000];
    double w2[3000][3000];
    r[s+1][n] = u2[s+1][n]*w2[s+1][n]*w2[s+1][n] - b*w2[s+1][n];
    return 0;
}

発生している問題・エラーメッセージ

gccは通るのですが、プログラムを実行すると強制終了してしまいます。

Answer 1

突っ込みどころはたくさんありますね。
ざっと気づいたところ。

    u2[s][n]=1.0;
    w2[s][n]=1.0;

sとnを初期化していないので、どこが初期化されるかわかりません。

    for(n=0;n<=3000;n++){
        for(s=0;s<=3000;s++){
            for(n=0;n<=3000;n++){

nのループが2回入れ子になってる。
ループ回数が0～3000になってる。
配列の要素は0～2999です。

        for(s=0;s<=3000;s++){
            for(n=0;n<=3000;n++){
                u2[s+1][n] = (mp*u2[s][n] + delta*(2*mp*u2[s][n]*w2[s+1][n]+a) )/1 + delta*(w2[s+1][n]*w2[s+1][n]+a+1);
                w2[s+1][n] = (mq*w2[s][n] + delta*(mp*u2[s][n]*(mq*w2[s][n]*w2[s][n]+1)+b) )/1 + delta*(w2[s+1][n]*w2[s+1][n]+a+1);
            }

u2[s+1][n]では、最大でu2[3001][3000]を参照することになり、確保していない領域を参照することになります。

    double k0[2][2], k1[2][2], k2[2][2], k3[2][2];//計算した値を入れる配列

コード内で使ってる要素は[0][0]か[1][1]なので、二次元配列にする意味がない。

コード全体を通して、[3000][3000]の配列。

初期化せずに使っている。
またサイズが大きすぎてスタックオーバーしているはず。
グローバル変数にするかstatic変数にしてください。

Comments

[S.Abe]

細かいところまでご指摘ありがとうございました。何とか修正できましたのでベストアンサーにさせていただきます。

Answer 2

関数内で3000x3000なdoubleの配列が３つも宣言されてますのでスタックオーバーフロー起こしてたりしませんか？
3000x3000って900万になりますけどスタック領域そんなに無いですよね？
mallocで確保するとかしなくて大丈夫です？

Answer 3

まず問題になりそうなのはここですかね。

int main(void) {
  FILE *output;
  output = fopen("output.dat","w");

  double t = 0.0;
  double dt = 0.01;
  double delta = 0.1;
  double r[3000][3000];
  double u = 1.0, w = 1.0;//初期値
  double a=0.03, b=0.10;//論文での値
  double k0[2][2], k1[2][2], k2[2][2], k3[2][2];//計算した値を入れる配列

  //離散化を構成する
  int i;
  int n,s,p;//mpとmqの添え字
  double mp,mq;
  double u2[3000][3000];
  double w2[3000][3000];
  //double u2[0][0]=0.0;//初期化
  //double w2[0][0]=0.0;//初期化
  u2[s][n]=1.0;　// ①
  w2[s][n]=1.0;

  for(n=0;n<=3000;n++){
    for(s=0;s<=3000;s++){
      for(n=0;n<=3000;n++){ // ②

まず①ですが、u2[s][n] にアクセスしようとしていますが、s, n が 初期化されていません。
デバッガを介した場合は 0 クリアされていることがありますが、デバッガを介さないと s, n に何が入っているか分かったものではありませんから、変な場所にアクセスして segmentation fault で落ちますな。

そして②ですが、一番大きいループで n を回しているのに、その中でまた n を回すって、何がやりたいんですか？