试图模拟双摆的分段错误

Question

我正在尝试编写一个程序来模拟双摆，就像this program。当我运行这段代码时，我遇到了一个段错误，但我找不到导致它的原因。

#include <iostream>
#include <cmath>
#include <vector>

#define PI 3.14159265 // Pi
#define G 9.8 // Acceleration due to gravity
#define L1 1.0 // Pendulum 1 Length ( meters )
#define M1 1.0 // Pendulum 1 Mass ( kilograms )
#define L2 1.0 // Pendulum 2 Length ( meters )
#define M2 1.0 // Pendulum 2 Mass ( kilograms )
#define N 4 // Number of equations. Don't change

using namespace std;

void calculateDerivatives ( vector<float> &dydx , vector<float> &yIn )
{
   // Calculate derivatives and put them in dydx vector
   dydx.clear();
   float temp1 , temp2 , delta;

   delta = yIn[2] - yIn[0];

   // Calculate sins and cosines to lower computation cost
   float sD = sin ( delta );
   float cD = cos ( delta );
   float s2D = sin ( 2 * delta );

   dydx.push_back ( yIn[1] );

   temp1 = ( M2 * L1 * pow ( yIn[1] , 2 ) * (1 / 2) * s2D );
   temp1 += ( M2 * G * sin ( yIn[2] ) * cD );
   temp1 += ( M2 * L2 * pow ( yIn[3] , 2 ) * sD );
   temp1 -= ( ( M1 + M2 ) * G * sin ( yIn[0] ) ) / ( ( M1 + M2 ) * L1 - M2 *L1 * pow ( cD , 2 ) );
   dydx.push_back ( temp1 );

   dydx.push_back ( yIn[3] );

   temp2 = ( -M2 * L2 * pow ( yIn[3] , 2 ) * ( 1 / 2 ) * s2D );
   temp2 += ( M1 + M2 ) * G * sin ( yIn[0] ) * cD;
   temp2 -= ( M1 + M2 ) * L1 * pow ( yIn[1] , 2) * sD;
   temp2 -= ( ( M1 + M2 ) * G * sin ( yIn[2] ) ) / ( ( L2 / L1 ) * delta );
   dydx.push_back ( temp2 );
}

void doRungeKutta ( const float x , const float h , vector<float> &yIn , vector<float> &yOut )
{
   // Use the Runge Kutta method to solve the DE
   vector<float> dydx;
   vector<float> dydxt;
   vector<float> yt;
   vector< vector<float> > kuttas;

   // Step 1
   calculateDerivatives ( dydx , yIn );
   vector<float> temp;
   for ( int i = 0; i < N; i++ )
   {
      temp.push_back( h * dydx[i] );
      yt.push_back ( yIn[i] + ( 0.5 ) * temp[i] );
   }
   kuttas.push_back ( temp );

   // Step 2
   calculateDerivatives ( dydxt , yt );
   for ( int i = 0; i < N; i++ )
   {
      temp[i] = ( h * dydxt[i] );
      yt[i] = ( yIn[i] + ( 0.5 ) * temp[i] );
   }
   kuttas.push_back ( temp );

   // Step 3
   calculateDerivatives ( dydxt , yt );
   for ( int i = 0; i < N; i++ )
   {
      temp[i] = ( h * dydxt[i] );
      yt[i] = ( yIn[i] + temp[i] );
   }
   kuttas.push_back ( temp );

   // Step 4
   calculateDerivatives ( dydxt , yt );
   for ( int i = 0; i < N; i++ )
   {
      temp[i] = ( h * dydxt[i] );
      yOut.push_back ( yIn[i] + kuttas[0][i] / 6.0 + kuttas[1][i] / 3.0 + kuttas[2][i] / 3.0 + kuttas[3][i] / 6.0 );
   }

}

int main( int argc , char *argv[] )
{
   int steps;
   float time_min , time_max , sTheta1 , sTheta2 , sAngVel1 , sAngVel2;
   float degRad = PI / 180;
   vector<float> yIn , yOut , times , theta1 , theta2 , angVel1 , angVel2;

   // Get command line inputs
   time_min = atof ( argv[1] );
   time_max = atof ( argv[2] );
   sTheta1 = atof ( argv[3] );
   sAngVel1 = atof ( argv[4] );
   sTheta2 = atof ( argv[5] );
   sAngVel2 = atof ( argv[6] );
   steps = atoi ( argv[7] );

   int h = ( time_max - time_min ) / ( steps - 1 );
   for ( int i = 0; i < steps; i++ )
   {
      times[i] = time_min + h * i;
   }

   theta1.push_back( sTheta1 * degRad );
   theta2.push_back( sTheta2 * degRad );
   angVel1.push_back( sAngVel1 );
   angVel2.push_back( sAngVel2 );

   cout << times[0] << " " << theta1[0] << " "
        << angVel1[0] << " " << theta2[0] << " " << angVel2[0];

   for ( int i = 0; i < ( steps - 1 ); i++ )
   {
      yIn.push_back( theta1[i] );
      yIn.push_back( angVel1[i] );
      yIn.push_back( theta2[i] );
      yIn.push_back( angVel2[i] );

      doRungeKutta ( times[i] , h , yIn , yOut );

      theta1.push_back( yOut[0] );
      angVel1.push_back( yOut[1] );
      theta2.push_back( yOut[2] );
      angVel2.push_back( yOut[3] );
      cout << times[i + 1] << " " << theta1[i + 1] << " "
           << angVel1[i + 1] << " " << theta2[i + 1] << " " << angVel2[i + 1];

   }

   return 0;
}

对不起，它太长了，我只是不知道是什么原因造成的，所以我无法将它缩短为有问题的部分。我一遍又一遍地检查代码，寻找我在声明数组或访问越界内容时可能犯的某种错误，但我找不到我的错误。

运行 gdb 我得到了这个：

$ g++ dp.cpp
$ ./a.out 0 10 90 0 90 0 1000
Segmentation fault (core dumped)
$ gdb a.out
GNU gdb (Ubuntu 8.1-0ubuntu3.2) 8.1.0.20180409-git
Copyright (C) 2018 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from a.out...(no debugging symbols found)...done.
(gdb) run 0 10 90 0 90 0 1000
Starting program: a.out 0 10 90 0 90 0 1000

Program received signal SIGSEGV, Segmentation fault.
0x0000000008001743 in doRungeKutta(float, float, std::vector<float, std::allocator<float> >&, std::vector<float, std::allocator<float> >&) ()
(gdb) where
#0  0x0000000008001743 in doRungeKutta(float, float, std::vector<float, std::allocator<float> >&, std::vector<float, std::allocator<float> >&) ()
#1  0x0000000008001c8a in main ()
(gdb)

这对我没有帮助。我尝试在 doRungeKutta 函数中搜索以找到问题，可惜我找不到。

编辑：使用向量而不是分配内存。但是，同样的 gdb 输出仍然会发生段错误。

Answer 1

这里有一个错误：

void doRungeKutta ( const float x , const float h , vector<float> &yIn , vector<float> &yOut )
{
   //...
   vector< vector<float> > kuttas;

   // Step 1
   calculateDerivatives ( dydx , yIn );
   //...
   kuttas.push_back ( temp );  // Here is the first one

   // Step 2
   calculateDerivatives ( dydxt , yt );
   //.. 
   kuttas.push_back ( temp );  // Here is the second one

   // Step 3
   calculateDerivatives ( dydxt , yt );
   //...   
   kuttas.push_back ( temp ); // Here is the third one

   // Step 4
   calculateDerivatives ( dydxt , yt );
   for ( int i = 0; i < N; i++ )
   {
      temp[i] = ( h * dydxt[i] );
      yOut.push_back ( yIn[i] + kuttas[0][i] / 6.0 + kuttas[1][i] / 3.0 + kuttas[2][i] 
                  / 3.0 + kuttas[3][i] / 6.0 );  // <-- Error -- there is no kutta[3]
   }
}

kuttas 向量的大小为 3，但您正在尝试访问 kutta[3]，这是一个越界访问。因此，内存覆盖和未定义的行为。

您还有其他错误，例如某些 vector 项目具有 -inf 值，但明显的错误是指出的错误。