1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90

#include <iostream>
#include <iomanip>
#include <valarray>
using namespace std;
// Function prototypes
void step( double dx, double &x, valarray<double> &Y );
valarray<double> F( double x, valarray<double> Y );
// Some global variables (hmm!  there are better ways) as kij is needed in F()
double k[2][2];
//======================================================================
int main()
{
int nstep;
double t, dt;
double x[2][2];
k[0][0] = 0.1; k[0][1] = 0.3; k[1][0] = 0.02; k[1][1] = 0.05; // no idea; I'm guessing
x[0][0] = 1.0; x[0][1] = 0.1; x[1][0] = 0.05; x[1][1] = 1.0 ; // ditto
// Set number of variables
int ndep = 4;
valarray<double> Y(ndep); // Take Y[0] ... Y[3] as X00, X01, X10, X11 (e.g.)
// Set initial values and global constants
t = 0; Y[0] = x[0][0]; Y[1] = x[0][1]; Y[2] = x[1][0]; Y[3] = x[1][1];
// Set step size and number of steps // probably going to need a smaller step size and more steps
dt = 0.1;
nstep = 10;
// Write header and initial conditions
#define SP << setw( 12 ) << fixed << setprecision( 4 ) << // Spacer: save some repetition when writing
cout << "t" SP "x00" SP "x01" SP "x10" SP "x11" << endl;
cout << t SP Y[0] SP Y[1] SP Y[2] SP Y[3] << endl;
// Solve the differential equation using nstep intervals
for ( int n = 1; n <= nstep; n++ )
{
step( dt, t, Y );
cout << t; // Rather crude output
for ( int i = 0; i < ndep; i++ ) cout SP Y[i]; // You probably want to dump to file instead
cout << endl;
}
}
//======================================================================
void step( double dx, double &x, valarray<double> &Y )
{
int ndep = Y.size();
valarray<double> dY1(ndep), dY2(ndep), dY3(ndep), dY4(ndep);
dY1 = F( x , Y ) * dx;
dY2 = F( x + 0.5 * dx, Y + 0.5 * dY1 ) * dx;
dY3 = F( x + 0.5 * dx, Y + 0.5 * dY2 ) * dx;
dY4 = F( x + dx, Y + dY3 ) * dx;
Y += ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
x += dx;
}
//======================================================================
valarray<double> F( double x, valarray<double> Y )
{
valarray<double> f( Y.size() ); // this will be returned as F in the output
double x00 = Y[0], x01 = Y[1], x10 = Y[2], x11 = Y[3]; // Makes it easier to code
double f00 = k[0][0] * x00 * ( 1 + x01 ) + k[0][1] * x00 * ( 1 + x00 ); // Not sure this is what you meant ????????
double f01 = k[0][1] * x01 * ( 1 + x00 ) + k[0][0] * x10 * ( 1 + x01 );
double f10 = k[1][0] * x10 * ( 1 + x11 ) + k[1][1] * x01 * ( 1 + x10 );
double f11 = k[1][1] * x11 * ( 1 + x10 ) + k[1][0] * x11 * ( 1 + x11 );
f[0] = f00;
f[1] = f01;
f[2] = f10;
f[3] = f11;
return f;
}
//======================================================================
 