Я написал программу для аппроксимации решений обыкновенных дифференциальных уравнений методом Адама.
Запуск программы с помощью 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);
Очевидно, что что-то не так с моей обработкой vector.push_back, но я не знаю, с чего начать.Я не могу вспомнить случай, когда изменение вектора является незаконным.
Вызовите дифференцировать (), чтобы начать.Математика делается в шаге ().Адаптивное продвижение времени через интервал с advance ().Проверьте выбранный шаг по времени с помощью checkdt (), прежде чем снова запустить step ().
Извините за огромный дамп кода.Я уверен, что многие улучшения могут быть сделаны исключительно с точки зрения 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";
}
Спасибо.