1

我正在尝试使用 Runge-Kutta 4 阶来解决双摆问题。当我尝试运行该程序时，它会崩溃，并且在循环一次后出现 c0000005 错误。我猜这可能与 System 函数中的长计算有关。来自更有经验的 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))
4

1 回答 1

2

该代码有一些突出的问题，您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。k正如乔纳森指出的那样，你的循环似乎表明你真的打算声明k[4][4]哪个可以解决这个问题。

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

于 2013-10-25T03:14:26.170 回答