Just transform your equations to a first order ODE and use a state type of length 2 N
. The first N entries now handle only the x[i]
while the second N entries refer to the velocities x'[i]
void ode( state_type const& x , state_type &dxdt , double t )
{
for( size_t i=0 ; i<N ; ++i )
{
double sum = 0.0;
// calculate sum
dxdt[i] = x[i+N];
dxdt[i+N] = K * sum;
}
}
A complete example might look like
size_t N = 512;
typedef std::vector< double > state_type;
state_type x( 2 * N );
// initialize x
double t_start = 0.0 , t_end = 10.0 , dt = 0.01;
odeint::integrate( ode , x , t_start , t_end , dt );