c - C、C0000005、および長い計算によるクラッシュエラーは初めてですか?

Question

二重振り子の問題を解決するために、ルンゲ クッタ 4 次を使用しようとしています。プログラムを実行しようとするとクラッシュし、ループを 1 回通過した後に c0000005 エラーが発生します。システム関数の長い計算と関係があるのではないかと推測しています。この問題の原因について、経験豊富な C ユーザーからのヒントはありますか?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//Attempt to solve double pendulum

//B,a to be used later. m=mass, l=length of rod, initial values for theta's and p's
double B,a,m=1,l=1, g=9.81,theta1i=1,theta2i=1,p1i=0,p2i=0;

static double System(int i,double theta1, double theta2, double p1, double p2);

int main()
{
FILE *f;
f = fopen("Question3.txt", "w");
int i,j;
double N=10000;
double h=1/N;

/*Runge Kutta*/
/*I need 3 sets of parameters k, thus each are defined as arrays of length 3.*/

double k[3][4],x[]={theta1i,theta2i,p1i,p2i}, t=0, T=200;
while(t<T)
{
    printf("%8.3lf\t%8.3lf\t%8.3lf\t%8.3lf\t%8.3lf\n",t, x[0],x[1],x[2],x[3]);
    fprintf(f,"%8.31f\t%8.31f\t%8.3lf\t%8.3lf\t%8.3lf\n",t, x[0],x[1],x[2],x[3]);

    for (j=0;j<4;j++)
    {
        for (i=0;i<4;i++)
        {
            if (j==0){k[i][0]=h*System(i,x[0], x[1], x[2], x[3]);}
            if (j==1){k[i][1]=h*System(i,x[0]+k[0][0]/2.0,x[1]+k[1][0]/2.0,x[2]+k[2][0]/2.0,x[3]+k[3][0]/2.0);}
            if (j==2){k[i][2]=h*System(i,x[0]+k[0][1]/2.0,x[1]+k[1][1]/2.0,x[2]+k[2][1]/2.0,x[3]+k[3][1]/2.0);}
            if (j==3){k[i][3]=h*System(i,x[0]+k[0][2],x[1]+k[1][2],x[2]+k[2][2],x[3]+k[3][2]);}
        }
    }

    for (i=0;i<4;i++)
    {
        x[i]=x[i]+(k[i][0]+2.0*k[i][1]+2.0*k[i][2]+k[i][3])/6.0;
    }
    t=t+h;
}
fclose(f);
return(0);
}

double System(int i,double theta1, double theta2, double p1, double p2)
{
if (i==0){return  (p1);}
if (i==1){return  (p2);}
if (i==2){return  ((-0.5*m*l*l)*((6/m*l*l)*(2*p1-3*cos(theta1-theta2)*p2)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))*(6/m*l*l)*(8*p2-3*cos(theta1-theta2)*p1)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))*sin(theta1-theta2)+3*g*sin(theta1)/l));}
if (i==3){return  ((-0.5*m*l*l)*(-(6/m*l*l)*(2*p1-3*cos(theta1-theta2)*p2)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))*(6/m*l*l)*(8*p2-3*cos(theta1-theta2)*p1)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))*sin(theta1-theta2)+g*sin(theta1)/l));}

}
//theta1dt=(6/m*l*l)*(2*p1-3*cos(theta1-theta2)*p2)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))
//theta2dt=(6/m*l*l)*(8*p2-3*cos(theta1-theta2)*p1)/(16-9*cos(theta1-theta2)*cos(theta1-theta2))