基于微分求解器的c++分割故障

我写了一个程序，用亚当法近似求解常微分方程。

使用gdb运行程序得到:

Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_PROTECTION_FAILURE at address: 0x00007fff5f3ffff8
0x0000000100003977 in std::vector<double, std::allocator<double> >::push_back (this=0x100005420, __x=@0x100005310) at stl_vector.h:604
604         this->_M_impl.construct(this->_M_impl._M_finish, __x);

显然我对向量的处理有问题。但我不知道从哪里开始找。我想不出哪一种情况下修改向量是非法的。

调用differentiate()开始。数学是按步骤()完成的。使用advance()在间隔上推进的自适应时间。在再次运行step()之前，用checkdt()检查所选的时间步长。

为大量代码转储感到抱歉。我确信，仅仅从c++的角度出发，不需要数学知识，也可以做出许多改进:

//============================================================================
// Description : An Adam's Method Ordinary Differential Equation Approximation
// Disclaimer  : Posted to StackOverflow for help on a segmentation fault
//============================================================================
#include <iostream> //IO
#include <vector> //std::vector
#include <cmath> //abs, pow, sqrt
#include <numeric> //accumulate
using namespace std;
/* Terminology:
 * f(t, y) = the ordinary differential equation that will be solved
 * y(t) = solution of f at point t.
 *  told = the previous point at which f was evaluated and solved
 *  tnow = the current point at which f is evaluated and solved
 *  tnew = the new (interesting) point at which f will be evaluated and solved
 *  
 *  Yold = the corrected solution of the differential equation at told
 *  Ynow = the corrected solution of the differential equation at tnow
 *  Ynew = the corrected solution of the differential equation at tnew
 *  
 *  fold = the value of the given differential equation at told
         = f(told, Yold)
 *  fnow = the value of the given differential equation at tnow
         = f(tnow, Ynow)
 *  fnew = the value of the given differential equation at tnew
         = f(tnew, Ynew)
 *
 *  Pnew = prediction for the value of Ynew
 *  dt = abs(tnew - tnow)
 *  dtold = abs(tnow - told)
 */
//Information storage
struct simTime {
    double told;
    double tnow;
    double tnew;
    double dt;
    double dtold;
    double tol;
    double agrow;
    double ashrink;
    double dtmin;
    double dtmax;
    double endTime;
    double fold;
    double fnow;
    double fnew;
    double Yold;
    double Ynow;
    double Ynew;
    double Pnew;
    int stepsSinceRejection;
    int stepsRejected;
    int stepsAccepted;
} test;
//Define global variables
std::vector<double> errorIndicators(0);
std::vector<double> solutions(0);
std::vector<double> differencesDDY(0);
std::vector<double> differencesDDYSquared(0);
std::vector<double> timesTNew(0);
//Function declarations
void step(double fnow, double fold, double Ynow, double tnew, double tnow,
        double dtold, double dt, double(*f)(double t, double y), double told);
void advance();
void shiftvariables();
void printvector(std::vector<double>& vec);
void differentiate(double(*f)(double t, double y), double y0, double a,
        double b);
double f(double t, double y);
void checkdt();
int main() {
    differentiate(f, 0, 1, 5);
    cout << "Time values:" << endl;
    printvector(timesTNew);
    cout << "Solutions:" << endl;
    printvector(solutions);
    cout << "Differences between Prediction and Solution:" << endl;
    printvector(differencesDDY);
    return 0;
}
//Shift back all the variables to make way for the new values
void shiftvariables() {
    test.tnow = test.tnew;
    test.dtold = test.dt;
    test.Yold = test.Ynow;
    test.Ynow = test.Ynew;
    test.fold = test.fnow;
    test.fnow = test.fnew;
    advance();
}
//Ordinary differential equation to be solved
double f(double t, double y) {
    return pow(t, 2);
}
//Calculate the predicted and corrected solution at a chosen tnew
void step(double fnow, double fold, double Ynow, double tnew, double tnow,
        double dtold, double dt, double(*f)(double t, double y), double told) {
    //The calculation for Ynew requires integration. I first thought I would need to
    //  use my project 1 code to calculate the integration, but now I see in class we
    //  solved it analytically such that integration is not required:
    //Linear prediction of Ynew using Ynow and fnow
    double Pnew = Ynow + (dt * fnow) + (dt * dt / (2 * dtold)) * (fnow - fold);
    test.Pnew = Pnew;
    //Predict the value of f at tnew using Pnew
    double fnew = f(tnew, Pnew);
    test.fnew = fnew;
    //Calculate the corrected solution at tnew
    double interpolationFactor = fnew - (fnow + dt * (fnow - fold) / dtold);
    double integration = (dt / 6) * (2 * dt + 3 * dtold) / (dt + dtold);
    double Ynew = Pnew + interpolationFactor * integration;
    test.Ynew = Ynew;
    //Update the variables for the next round
    shiftvariables();
}
//Check the previous solution and choose a new dt to continue evaluation
void advance() {
    //The error indicator is the l2-norm of the prediction minus the correction
    double err_ind = sqrt(
            std::accumulate(differencesDDYSquared.begin(),
                    differencesDDYSquared.end(), 0));
    errorIndicators.push_back(err_ind);
    // Case where I reject the step and retry
    if (err_ind > test.tol && test.dt > test.dtmin) {
        ++test.stepsRejected;
        test.stepsSinceRejection = 0;
        test.dt = test.dt * 0.5;
        test.tnew = test.tnow + test.dt;
        checkdt();
    }
    // Cases where I accept the step and move forward
    else {
        ++test.stepsAccepted;
        ++test.stepsSinceRejection;
        solutions.push_back(test.Ynew);
        differencesDDY.push_back(abs(test.Pnew - test.Ynew));
        differencesDDYSquared.push_back(pow((test.Pnew - test.Ynew), 2));
        //Decrease dt
        if (err_ind >= 0.75 * test.tol) {
            test.dtold = test.dt;
            test.dt = (test.dt * test.ashrink);
            test.tnew = test.tnow + test.dt;
            checkdt();
        }
        //Increase dt
        else if (err_ind <= test.tol / 4) {
            if ((test.stepsRejected != 0) && (test.stepsSinceRejection >= 2)) {
                test.dt = (test.dt * test.agrow);
                test.tnew = test.tnow + test.dt;
                checkdt();
            } else if (test.stepsRejected == 0) {
                test.dt = (test.dt * test.agrow);
                test.tnew = test.tnow + test.dt;
                checkdt();
            }
        }
    }
}
//Check that the dt chosen by advance is acceptable
void checkdt() {
    if ((test.tnew < test.endTime) && (test.endTime - test.tnew < test.dtmin)) {
        cout << "Reached endTime." << endl;
    } else if (test.dt < test.dtmin) {
        test.dt = test.dtmin;
        test.tnew = test.tnow + test.dt;
        timesTNew.push_back(test.tnew);
        step(test.fnow, test.fold, test.Ynow, test.tnew, test.tnow, test.dtold,
                test.dt, f, test.told);
    } else if (test.dt > test.dtmax) {
        test.dt = test.dtmax;
        test.tnew = test.tnow + test.dt;
        timesTNew.push_back(test.tnew);
        step(test.fnow, test.fold, test.Ynow, test.tnew, test.tnow, test.dtold,
                test.dt, f, test.told);
    } else if ((test.tnew + test.dt) > test.endTime) {
        test.dt = test.endTime - test.tnew;
        test.tnew = test.tnow + test.dt;
        checkdt();
    } else if (((test.tnew + test.dt) < test.endTime)
            && ((test.tnew + 2 * test.dt) > test.endTime)) {
        test.dt = (test.endTime - test.tnew) / 2;
        test.tnew = test.tnow + test.dt;
        checkdt();
    }
    //If none of the above are satisfied, then the chosen dt
    //  is ok and proceed with it
    else {
        timesTNew.push_back(test.tnew);
        step(test.fnow, test.fold, test.Ynow, test.tnew, test.tnow, test.dtold,
                test.dt, f, test.told);
    }
}
//meta function to solve a differential equation, called only once
void differentiate(double(*f)(double t, double y), double y0, double a,
        double b) {
    //Set the starting conditions for the solving of the differential equation
    test.fnow = f(a, y0);
    test.endTime = b;
    test.Ynow = y0;
    //solutions.push_back(y0);
    timesTNew.push_back(a);
    //Set the constants
    test.ashrink = 0.8;
    test.agrow = 1.25;
    test.dtmin = 0.05;
    test.dtmax = 0.5;
    test.tol = 0.1;
    //Set fold = fnow for the first step
    test.fold = test.fnow;
    test.tnow = a;
    test.told = a - test.dtmin;
    test.dtold = abs(test.tnow - test.told);
    //Create the first prediction, which will then lead to correcting it with step
    advance();
}
// Takes a vector as its only parameters and prints it to stdout
void printvector(std::vector<double>& vec) {
    for (vector<double>::iterator it = vec.begin(); it != vec.end(); ++it) {
        cout << *it << ", ";
    }
    cout << "n";
}

谢谢。