4th order Runge-Kutta method of vectors

[Kloppite]

Write your question here.
Hi, guys, kindly chip in with advice regarding how I should go about this program. I have intermediate knowledge of c++ (self0taught) but my current version of this program while effective appears to take a while to run for some functions.

Basically, i am trying to solve 10 coupled differential equations using the 4th order RK method. The first portion of my code opens data files for the initial conditions of the 10 variables and some parameters that I am using.

The next portion of my code should be for the iteration but I am sort of stuck here since initially, I have written the code for four of the equations in this form


what I find though is that the function takes too many arguments which make the code looks all jammed up. Will appreciate any advice on implementation and functionality as

[lastchance]

Use std::valarray for your dependent variable Y and for its derivative F. Please!

Then you can write and solve a system that looks almost indistinguishable from the scalar version:
dY/dx = F(x,Y)

Write a function of the form
valarray<double> F( double x, valarray<double> Y )
to return your derivative function. This is the only part which depends on the equation being solved.

Then you won't have to write 4 lines (or worse) for each step.

Here's a version of one step. All you have to do is write the function F for your floating device (or whatever it is).

void step( int ndep, double dx, double &x, valarray<double> &Y, valarray<double> (*F)( double, valarray<double> ) )
{
   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;
}

[Kloppite]

Thanks for the input, but this is a coupled system and each equation is linked somehow. I cannot iterate over one without considering the contributions from the other hence the whole cluster of a code I had up there. Do you suggest an alternative?

[lastchance]

This is a standard way of solving coupled equations.

At each step Y[0], Y[1], Y[2], Y[3] are your current dependent variables and are the components of a 4-element valarray.

Simply write the derivative function components F[0], F[1], F[2], F[3] in terms of these and return it as a valarray F.

I'm assuming from your example that ndep, the number of dependent variables, which are presumably positions or angles, is 4 ... but the step routine would work for anything. Just code the appropriate derivative function.

[Kloppite]

Where I seem to not understand is the void function you have above, as it looks to me like I can only get single element iterations over each one. Will there be a for loop somewhere before the void function or will the valarray take care of that also?

Thanks.

[lastchance]

That void function is just doing one complete (time)step or whatever. Can you see that it has precisely the same 4 steps as you have? And the same weighted average at the end? It's just that the substeps for the different variables are associated with the different components of the valarray ... and the use of a valarray allows you to do vector addition just like you would in maths as a single equation, not having to have a separate line for each component.


I tried to write you an example of how you would use this technique to solve a slightly simpler set of coupled equations ... just 2 of them (you have 4, I think).

Consider a simpler system with just TWO (ndep=2) coupled equations:

du                        dv
-- = x + u - 3v           -- = 2x - 6u + v
dx                        dx

Now you could write these as the single VECTOR equation

d |u|   | x + u - 3v|   
--| | = |           |   
dx|v|   |2x - 6u + v|   

But, if you write Y0=u and Y1=v then these are

d |Y0|   | x +  Y0 - 3Y1|
--|  | = |              |
dx|Y1|   |2x - 6Y0 +  Y1|

or, simply

d |Y0|   |F0|               dY                           
--|  | = |  |          i.e. -- = F(x,Y)
dx|Y1|   |F1|               dx

where

F0 =  x + Y0 - 3Y1      i.e. a function of x and Y
F1 = 2x - 6Y0 + Y1      i.e. another function of x and Y

Can you see that we simply have to write the components of the derivative function (F0,F1) in terms of x and the original dependent variables (Y0,Y1), provided we make the association of Y0 with u and Y1 with v.


Your problem is a little bigger ... but only because it has 4 components instead of 2.

[Kloppite]

Thanks for the explanations above.

I have declared a Valarray function like this

valarray<double> F(double x, valarray<double> Y){
return F(x, Y);
}

How do i pass the vector values into it, say i wanted to pass F0 from the equation above, since it is a function of Y[1], Y[2] and x, do i explicitly write is as such?
e.g

F(x[i], Y[1][0], Y[2][0])

[lastchance]

Here is an example of how you would solve my 2-component example

du                        dv
-- = x + u - 3v           -- = 2x - 6u + v
dx                        dx

with a somewhat arbitrary initial condition u=2.5, v= -1.2 at x = 0.

Note that I have simplified the step function slightly (the original came from a bigger project) by dropping ndep and the derivative function, provided that I define F() in the rest of the program. ndep can get picked up automatically from the number of components of Y.

See how the derivative function is encoded. (Note that it comes back as the name of the function, upper-case F. You actually need a different variable within the function, so I have denoted this different variable lower-case f for convenience.)

#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 );

//======================================================================


int main()
{
   int nstep;
   double x, dx;

   // Set number of variables
   int ndep = 2;
   valarray<double> Y(ndep);

   // Set initial values
   x = 0;   Y[0] = 2.5;   Y[1] = -1.2;

   // Set step size and number of steps
   dx = 0.1;
   nstep = 10;


   // Write header and initial conditions
   #define SP << setw( 12 ) << fixed << setprecision( 4 ) <<    // save some repetition when writing
   cout << 'x' SP 'u'  SP 'v'  << endl;
   cout <<  x  SP Y[0] SP Y[1] << endl;

   // Solve the differential equation using nstep intervals
   for ( int n = 1; n <= nstep; n++ ) 
   {
      step( dx, x, Y );                          // main Runge-Kutta step
      cout << x SP Y[0] SP Y[1] << endl;         // output
   }
}

//======================================================================

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

   f[0] =       x +     Y[0] - 3 * Y[1];         // ****************************************
   f[1] = 2.0 * x - 6 * Y[0] +     Y[1];         // just encode derivative components here *
                                                 // ****************************************
   return f;
}

//====================================================================== 

x           u           v
0.0000      2.5000     -1.2000
0.1000      3.4289     -3.1450
0.2000      5.2207     -6.1064
0.3000      8.4005    -10.8831
0.4000     13.8737    -18.7967
0.5000     23.1840    -32.0595
0.6000     38.9481    -54.3937
0.7000     65.5927    -92.0748
0.8000    110.5991   -155.6944
0.9000    186.6061   -263.1342
1.0000    314.9611   -444.5902

Your derivative function will be more complex to code, obviously, but you do so in the same place.


How complicated would it be for you to actually write down your original differential equations here? Also, can you state what the physical problem is and what your dependent variables are. (Guessing at positions, angles, maybe velocities?)

[Kloppite]

Thanks for the explanations. It does seem a little clearer now. It's actually a set of compartmental flow equations so the dependent variables are concentrations in time. Say, we have flow between 2 compartments/positions 1&2, of two different quantities A&B, my differential equations represents


$dx[i][j]/dt = k[i][j]x[i][j](1+x[i][j']) + k[i][j']x[j][i](1+x[i][j])$

Here;
i's are the quantities A,B.
j= initial compartment/position.
j' = final compartment/position.
k = velocity depending on the i/j/j' although it is a constant.
x[i][j][n] = concentration of quantity (i), in compartment(j), step(n).


From this example 1 have four differential equation since I have to determine x[1][1][n],x[2][1][n],x[1][2][n],x[2][2][n].

My program (from which the portion in the OP was taken) actually has a lot more (i's and j's, and j''s) but I have written a mess of a code which works well but takes a while to run as it is a sea of conditional statements and loops, hence my question.

Thanks.

[lastchance]

OK, I'm not sure that I'm interpreting your equations correctly, but I hope that you can see roughly what I'm doing. I'm not convinced your equations would look exactly like that as concentrations would blow up exponentially without some minus signs somewhere.

Note then for the following:
- I'm not quite sure about your equations;
- I've mapped Y[0] ... Y[3] to X00, X01, X10, X11
- I've changed the independent variable from x to t (for time) in main, because (a) it's more physical and (b) otherwise there would be confusion with xij
- you should not store things at all time levels [n]; you should write concentrations to file; at the moment the output is crude and to screen, but you will need to output it in whatever form you require.

So this is roughly what you seem to need (but please fix those equations):

#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;
}

//====================================================================== 

t         x00         x01         x10         x11
0.0000      1.0000      0.1000      0.0500      1.0000
0.1000      1.0748      0.1069      0.0507      1.0093
0.2000      1.1580      0.1145      0.0515      1.0187
0.3000      1.2511      0.1229      0.0524      1.0282
0.4000      1.3557      0.1323      0.0533      1.0379
0.5000      1.4742      0.1429      0.0542      1.0476
0.6000      1.6094      0.1548      0.0552      1.0575
0.7000      1.7647      0.1685      0.0563      1.0675
0.8000      1.9452      0.1842      0.0575      1.0776
0.9000      2.1570      0.2026      0.0587      1.0878
1.0000      2.4090      0.2243      0.0601      1.0982

[Kloppite]

Hi, Thanks a lot.

I have not written out my equation in its whole complexity actually that's why it appears to be a little crude and basic. I do have conditions to store only positive values and expected errors are also considered. I have also been writing/reading final/initial values of variables and parameters to and from files, as well as a convergence function.

My initial program with all the conditions and function is about 810 lines - making these adjuments will maybe make it a lot more compact and effective.

Once again, thank you for your time. I will update the topic as I proceed along but you have been of great help in understanding something I reckon I have a great need of.

[Kloppite]

Just another point, I never knew I could have arrays as global variables. I have about 20 such parameters in my cluster and sometimes I have written functions with up to 15 4-D arrays.

[48T7M4Gy]

.

[Kloppite]

Hey, Guys how do you reckon I need advice on a number of things.

First, I have a number of input files that I open and read values from within the main function, however, I will like to open and read these files from a different function, which when called from main returns all the parameter and variables(stored initially in a multidimensional array).

Secondly, in solving the numerical equations that I have (14 in total), I do have to find the convergence of some of them, how reasonable is it to write another function for convergence, which again takes in the prior values of my variable in order to give a new converged value.

Please keep in mind I have written this code before with the convergence included but not within a function, but with a whole chunk of messy equations just as I had on there initially.

Thank you.

@kemort @lastchance

[48T7M4Gy]

.

[Kloppite]

Thanks for the prompt response. So, my dimensions are
i=Item(I have 5 items)
j=location-in(I have 4 location-in)
k=location-to(I have 4 of these also)
l=process-number(three of these)

Now these indices are for parameters and they are constants, about 4 of them which are unique to specific items in specific locations.

Hence, I can say I have T[5][4][4][3] as a parameter for example.

The dependent variables themselves are in the form

X[il[j](t+nh)

Again with 5 i's and 4 J's.

[lastchance]

Your description is just too cryptic: objects with 4 dimensions are few and far between. I suspect now that your DEPENDENT variables are:
- concentration
- species
In other words, X[species i][location j] - that is TWO dimensions.

What you are describing above seems to be the TERMS which COUPLE those equations, which could actually be read, stored and used as a ONE-dimensional array of structs; e.g.
struct process{ int species, in, out; double value; }. When you read them, just assign them sequentially to a vector of structs using push_back(). When you use them in the derivative function just go through them one by one.

If , as I suspect, t is time, you should not be storing that - you should dump your X to file at identified intervals.

It is unclear what you mean by "convergence". This version of RK is explicit: there is no in-step iteration. If you meant convergence to final steady state then just use the normed distance between successive values of your dependent-variable array.
I think you need to write down your original equations.

[48T7M4Gy]

.

[Kloppite]

Thanks Kemort and lastchance for the ideas.

I will try to explain what I am trying to do as much as I can.

I have a set of variables which change in time governed by a set of differential equations. These variables I will call X[I][j] are are set of concentration, x, of item, I, in position j. I have 5 items and 4 positions.
Hence I am looking for 20 different concentrations in time: X[5][4]

The differential equations are of this form:
dx[I][j]/dt =-k[i][j][l][n]x[I][j](1+x[i][l]) + k[i][l][j][n]x[I][l](1+x[i][j]) and some extra paramaters multiplied by variables in this form.

The parameter k[i][j][l][n] here is a rate constant and it describes the rate of flow of item I, from position j, to l, in kinetic process (reaction order (n = 0 to 2)). Likewise for other parameters which have this form.

Hitherto, I have written one mess of a code for the case where I=1, that is four variables x[1][4] with my basic c++ knowledge, which gives accurate results, but is a mess of a code with lots of conditionals and takes a while to run even if it gives accuracy when compared with theoretical values.

In my code I had the following syntax:

1. Opened the files (10) in total
2. Made the parameter into 4D arrays as described above
3. Added an extra 'step' dimension to the variable array so that I have x[1][j=position][n=iteration-step]
4. Calculated the 4-Rk steps for each of these variables making 16 equations(k1..K4) in total.
5. Solved for x[1][j][n+1]
6. Did convergence by reducing h within a while loop which again had all sixteen equations :((
7. Solve for converged values of each which is then assigned to x[1][j][n+1]
8. Wrote the output to files
9. Closed the files.
I know this is the entirely wrong way of doing it but I am very stuck right now. I have been at it for the best part of three months. Any advice is appreciated.

[48T7M4Gy]

.