新的C，C0000005和碰撞错误，因为计算时间长？

Question

我试着用龙格-库塔第四阶来解决一个双摆问题。当我试图运行这个程序时，它会崩溃，在循环一次之后，我会得到一个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))

Answer 1

代码中有一些突出的问题，您在几个地方访问k是越界的，例如：

 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);}
             ^^^

k被声明为k[3][4]，因此3是k界限之外的一个，因为起诉开始于0。正如Jonathan所指出的，您的循环似乎表明您确实打算将k声明为k[4][4]，这将解决此问题。

在System中，如果i不是0 to 3，那么您将在没有返回的情况下退出函数，这是值返回函数的未定义行为。正如乔纳森已经提到的，System可能不是一个很好的名字，因为它非常接近标准的system()函数。