c++ - 如何以正确的方式解方程并保存？

Question

我是 C++ 的新手。所以实际上我会尝试找出每个时间步 (dt=0.00001) 的总力、速度、位置和总时间，并通过求解简单方程保存 myfile1、myfile2、myfile3、myfile4。我面临着找出确切值(value)的问题。当我运行程序时，发现段错误。那么我该如何解决这个问题呢？我附上了下面给出的程序:任何人都可以帮助我吗？

//C++ programming for selection random number
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <ctime>
#include <sstream>
#include <string>
#include <iomanip>
#include <cmath>


//%stifness
double kp=1000;
//%gravitational
double g=-9.81;


    double x[10000], y[10000], z[10000];
    double Fx[10000],Fy[10000],  Fz[10000];
    double vx[10000],vy[10000], vz[10000] ;
    double ax[10000],ay[10000], az[10000] ;
    double force[10000];
    const double PI=4.0*atan(1.0); 

    using namespace std;
int main()
{
    srand((unsigned) time(NULL));
//open output file
    ofstream myfile1;
    myfile1.open("1.dat");
    ofstream myfile2;
    myfile2.open("2.dat");
    ofstream myfile3;
    myfile3.open("3.dat");
    ofstream myfile4;
    myfile4.open("4.dat");


// %particle parameter

double  dt=1e-5;
double Ntimestep=1000;

//particle properties
    double R=0.003;
    double Dens=2600;
        int npart=1;

//Particle Creation 
        double Mass=(Dens*4*PI*(pow(R,3))/3);    //m=(4/3)* Dens*Pi*R^3
        cout<<"Mass="<< Mass<<endl;

//initial position of the particle

        x[0]=0;
            y[0] =0.01;
            z[0] =0;

// movement of the particle
// Particle initial velocity and forces
     Fx[0]=0;
     Fy[0]=0; 
     Fz[0]=0;
     vx[0]=0;
     vy[0]=-0.5;
     vz[0]=0;
     force[0]=0;

// Relation between dashport coefficient and coefficient of restitution
     double exp=0.9;
     double lnexp=0.10536;
     double Eta= ((2*sqrt(Mass*kp)*lnexp)/(sqrt(pow(PI,2)+pow(lnexp,2))));

//Time step
    int t=0;
for (int i=0;i<Ntimestep;i++)
{

// calculate particle to wall contact force

       if (y[i]<R)
            {
                Fy[i]=(-kp*(y[i]-R))-Eta*vy[i];
            }

// Calculate initial acceleration

        ay[i]=(Fy[i]/Mass)+g;
//force[i]+=force[1];
        force[i+1]=ay[i]*Mass;
        //cout<<"Total_Force="<<force<<endl;

// update valocity and displacement/location

    vy[i+1]=(vy[i]+(ay[i]*dt));
    y[i+1]=y[i]+vy[i]*dt+0.5*ay[i]*pow(dt,2);
    Fy[i]=0;

    t=t+1;

     double time=t*dt;



//...............output/save file..............................

cout<<"Total force="<<force<<endl;
myfile1<<"Total force="<<force<<endl;

cout<<"velocity="<<vy<<endl;
myfile2<<"velocity="<<vy<<endl;

cout<<"location="<<y<<endl;
myfile3<<"location="<<y<<endl;

cout<<"Total time"<<time<<endl;
myfile4<<"Total time="<<time<<endl;

}


//system ("PAUSE");
    cin.ignore();
    cin.get();
    myfile1.close();
    myfile2.close();
    myfile3.close();
    myfile4.close();



    return 0;
}

Answer 1

I suppose that your original value for Ntimestep was 10000 (the dimension of your C-style arrays x , y , z , Fx , Fy , Fz , vx , vy , vz and force ), not the actual 1000 .

对于实际的 Ntimestep = 1000 ，我没有看到出现段错误的情况；但是对于 Ntimestep = 10000 ，当 i 循环中的 for 为 9999(最后一次迭代)时，在以下说明中

force[i+1]=ay[i]*Mass;
vy[i+1]=(vy[i]+(ay[i]*dt));
y[i+1]=y[i]+vy[i]*dt+0.5*ay[i]*pow(dt,2); 

你在 force 、 vy 和 y 的位置 10000 中写入。这很糟糕(你是 C++ 的学习者，但我想你知道你可以读/写一个从 X x[n] 到 0 的 n-1 数组，并且访问在 x[n] 中是越界访问)。

一些建议，排名不分先后

1) 注意缩进；你的代码很难阅读

2) 尽可能避免使用全局变量

3) 尽可能使用常量(如果您使用 C++11/C++14，则尽可能使用 constexpr)；全局常量很好，您可以用它们定义数组大小；例如，您可以定义一个全局常量

 const int  numSteps = 10000;

和数组定义(我建议在 main() 本地)可以是

 double x[numSteps+1], y[numSteps+1], z[numSteps+1];
 double Fx[numSteps+1],Fy[numSteps+1],  Fz[numSteps+1];
 double vx[numSteps+1],vy[numSteps+1], vz[numSteps+1] ;
 double ax[numSteps+1],ay[numSteps+1], az[numSteps+1] ;
 double force[numSteps+1];

(+1给数组另一个位置，避免上面的问题)和for循环可以

 for (int i=0;i<numSteps;i++)

4) 考虑假设使用 std::vector<double>(带有 resize(numSteps+1))而不是 C 风格的数组；使用 std::vector ，使用 at() 而不是 [i] ，你会有一个异常，你会立即发现越界错误

4 之二)如果您使用的是 C++11/C++14，std::array<double, numSteps+1> 可能是比 std::vector<double> 更好的解决方案

p.s.:抱歉我的英语不好