CODEWITHSUNDEEP
vectordelaydifferential-equations
John Wiesel
John Wiesel

Reputation: 31

Solving DDE system

im trying to solve a differentiel delay equation system with c++. Im a newbie in terms of coding, so please if you have recommendations, tell me, I would like to improve my writing! What i want to do: initialize the history-array and then start to solve the differential equation by overwriting the history-array. But the problem is, I get the error message:

terminate called after throwing an instance of 'std::out_of_range' what(): vector::_M_range_check: __n (which is 9999) >= this->size() (which is 9999)

It seems that the history-arrays are out of range. I tried to put a std::cout in after the second if-condition to check if the code is going through the second for-loop, but he isn't. Since im learning c++ by doing right now, the problem isn't really clear to me. I hope someone sees the error. And dont hesitate to improve my code, I would really appreciate!

Thanks for your help!

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>

const double pi = 3.14159265358979323846;

//delay
int tau = 1;

//initial values
double x = 1.0;
double y = 1.0;
double t = 0.0;

//constants and parameters
double K = 0.25;
double lam = 0.5;
double omega = pi;

double dx, dy;

//step-size
double dt = pow(10.0, -4.0);

//number of steps
int Delta = static_cast<int>(tau/dt);

std::vector<double> hist_x((static_cast<int>(tau/dt) - 1), 0.0);
std::vector<double> hist_y((static_cast<int>(tau/dt) - 1), 0.0);

std::vector<double> t_val;
std::vector<double> x_val;
std::vector<double> y_val;

double euler(double f, double di, double time_step){
    f = f + time_step * di;
    return f;
    }

int main()
{
std::ofstream file_x;
std::ofstream file_y;
std::ofstream file_t;
file_x.open("x_val.txt");
file_y.open("y_val.txt");
file_t.open("t_val.txt");

for(int n = 0; n < 2; n++){
    if(n==0){
        for(int j; j < Delta; j++){
        dx = lam * x + omega * x;
        dy = lam * y - omega * x;

        x = euler(x, dx, dt);
        y = euler(y, dy, dt);
        t = t + dt;

        x_val.push_back(x);
        y_val.push_back(y);
        t_val.push_back(t);

        hist_x.at(j) = x;
        hist_y.at(j) = y;

        file_x<<x_val.at(j)<<std::endl;
        file_y<<y_val.at(j)<<std::endl;
        file_t<<t_val.at(j)<<std::endl;
        }
    }
    if(!(n==0)){
        for(int k = 0; k < Delta; k++){

        //f1(x,y)
        dx = lam * x + omega * x - K * ( x - hist_x.at(k) );

        //f2(x,y)
        dy = lam * y - omega * x - K * ( y - hist_y.at(k) );

        x = euler(x, dx, dt);
        y = euler(y, dy, dt);
        t = t + dt;

        x_val.push_back(x);
        y_val.push_back(y);
        t_val.push_back(t);

        hist_x.at(k) = x;
        hist_y.at(k) = y;

        file_x<<x_val.at(k + n * Delta)<<std::endl;
        file_y<<y_val.at(k + n * Delta)<<std::endl;
        file_t<<t_val.at(k + n * Delta)<<std::endl;
        }

    }


    }
file_x.close();
file_y.close();
file_t.close();
}

Upvotes: 0

Views: 79

Answers (1)

Armali
Armali

Reputation: 19375

        for(int j; j < Delta; j++){

You forgot to initialize j; you meant:

        for (int j = 0; j < Delta; j++)
        {

int Delta = static_cast<int>(tau/dt);

std::vector<double> hist_x((static_cast<int>(tau/dt) - 1), 0.0);
std::vector<double> hist_y((static_cast<int>(tau/dt) - 1), 0.0);

You index from 0 to Delta−1, this means the vectors need to have Delta elements, and you allocate one less; correct:

std::vector<double> hist_x(Delta, 0.0);
std::vector<double> hist_y(Delta, 0.0);

Upvotes: 2

Related Questions