C++ program has stopped working- Solving ordinary

Question

I'm writing a C++ program to find solutions for first order differential equations for a college assignment. The program starts up and then once I enter the number of iterations to do I get the error message "Euler's method.exe has stopped working". This is my code:

#include <functional>
#include <vector>

using namespace std;

 
double f_r(double x, double r) {  
  return r;
  } 

double f_s(double x, double s) {
  return -x/s;  
  }
 

double eulerstep(const function<double(double,double)>& f, double xsub0, double ysub0, double h) {
  double ysub1 = ysub0+ h * f(xsub0,ysub0);
  return ysub1;
  }


double euler(const function<double(double,double)>& f, double xsub0, double ysub0, double h, int n) {
  vector<double> xsub;
  vector<double> ysub;
  xsub[0] = xsub0;
  ysub[0] = ysub0;
  n = ysub.size();
  for (int i=1; i<n; i++){
  	xsub[i+1] = xsub[i] + h;
  	ysub[i+1] = ysub[i] + h * f(xsub[i],ysub[i]);
  	cout << xsub[i] << " , " << ysub[i] << endl;
  }
  return ysub[n];
  }
  
int main() {
  int nsteps = 0;
  cout << "Number of steps?" << endl;
  cin >> nsteps;
  double h = 1.0 / nsteps;

  double r = euler(f_r,0,1,h,nsteps);
  

  double s = euler(f_s,0,1,h,nsteps);
 
  return 0;
}

I'm fairly sure the problem lies in how I've defined my vectors but I'm new to using them so can't see where I've gone wrong. I'd be extremely grateful if anyone could point out the mistakes in my method

Thanks

Answer 1

vector<double> xsub;
vector<double> ysub;

You are instantiating a pair of vectors. They are initially empty.

xsub[0] = xsub0;
ysub[0] = ysub0;

You then proceed to assign values to the contents of the vector. This is where you crash, because the vectors are empty, and contain no values.

vector[x] references an existing element in a vector. The vector must already have at least x+1 elements in them, but none of your vectors have anything in them. They do not have element 0, nor element 1, nor any element.

Looks like your code expects each vector to contain n+1 elements, so you should explicitly invoke each vector's resize() method, accordingly, before you attempt to use each vector.

Answer 2

Here's an updated function. You have to set the sizes, not ask an empty vector for its size.

double euler(const function<double(double, double)>& f, double xsub0, double ysub0, double h, int n)
{
    vector<double> xsub;
    vector<double> ysub;
    xsub.resize(n+1); // so we can access [n], it must be size n+1
    ysub.resize(n+1);
    xsub[0] = xsub0;
    ysub[0] = ysub0;
    for (int i = 1; i<n; i++) {
        xsub[i + 1] = xsub[i] + h;
        ysub[i + 1] = ysub[i] + h * f(xsub[i], ysub[i]);
        cout << xsub[i] << " , " << ysub[i] << endl;
    }
    return ysub[n];
}