试图模拟双摆的分段故障

Question

我正在尝试制作一个模拟双摆的程序，很像这个程序。当我运行这个代码时，我得到了一个seg错误，但我找不到它的原因

#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函数中搜索以找到问题，可惜我找不到。

编辑：使用矢量而不是分配内存。然而，segfault仍然会在相同的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]，这是一种越界访问。因此，内存覆盖和未定义的行为。

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