2

I want to solve the same ODE but with several different initial conditions. I can solve with one fixed initial condition, for instance, in the example I have chosen x[0]=x[1]=x[2]=10.0. How could I create a while loop to solve the ODE for x[0]=10.0,9.0,....1.0? This is the code:

#include <vector>
#include<boost/numeric/odeint.hpp>

typedef std::vector<double> state_type;

const double sigma = 10.0;
const double R = 28.0;
const double b = 8. / 3.;

void lorenz(state_type &x, state_type &dxdt, double t)
{

    dxdt[0] = sigma * (x[1] - x[0]);
    dxdt[1] = R * x[0] - x[0] - x[0] * x[2];
    dxdt[2] = x[0] * x[1] - b * x[3];

}

int main()
{

    using namespace std;
    using namespace boost::numeric::odeint;
    state_type x(3);
    x[0] = x[1] = x[2] = 10.0;
    const double dt = 0.01;

    auto stepper = make_dense_output<runge_kutta_dopri5<state_type> >(1.0E-12,
                                                                      1.0E-12);
    double t = 0.0;
    stepper.initialize(x, t, dt);
    cout << t << " " << x[0] << endl;
    t = t + dt;
    while (t < 10)
    {
        if (t > stepper.current_time())
        {
            stepper.do_step(lorenz);
        }
        else
        {

            stepper.calc_state(t, x);
            cout << t << " " << x[0] << endl;
            t += dt;
        }

    }
    return 0;
}

Besides that, could anyone show me a minimum example or documentation where I can read about parallelize such ODE problem? I apologize for any inconvenient, I'm getting started in C++.

user4581301
  • 33,082
  • 7
  • 33
  • 54
Ninpou
  • 33
  • 3
  • What is your aim? Having an array or matrix of steppers that you then advance simultaneously? Or solving for the different initial conditions sequentially? Or making each integration a task which can be handed over to some thread pool, and organizing the logistics for the handling of the parallelism and the results? – Lutz Lehmann Dec 09 '19 at 20:22
  • Did you try to put the integration loop inside a loop over the `x[0]` values, with `stepper` as a block-local variable? Then this object is automatically deleted at the end of the loop and then constructed anew in the next loop. – Lutz Lehmann Dec 09 '19 at 20:26
  • I want to integrate the same ODE with different initial conditions and save the last value of x[0] for each initial condition into a matrix. I think what I need is making each integration a task which can be handed over to some thread pool and organizing the logistics for the handling of the parallelism and the results. – Ninpou Dec 09 '19 at 20:41
  • Since I am a novice in C++, I am trying to write the code using concurrent computation first (that is why I need the while loop for x[0]). However, I will need to parallelize this in the future because my real problem handles with ~10^6 values, that is, the dimension of my final matrix is ~10^3x10^3 – Ninpou Dec 09 '19 at 20:48

0 Answers0