4th order Runge-Kutta method of vectors

Pages: 123

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

Can you clarify that with regards to n? Are there SEVERAL terms on the rhs, one for each n, or what.

As @Kemort says, can you direct us to some weblink with the equations written out properly.

Kloppite (30)

@Kemort and @Lastchance

Here are some links for you sirs.

First one here is a prototype of my differential equations (not the exact one but close enough):

1. https://imgur.com/a/rpUou

Here, i is my present item(see below), j is my current compartment.
K[i][j][k'][l] = The rate of flow of item 'i' from compartment 'j', to compartment 'k' using process 'l'..Here, l determines A and B which can be an interger between 0 and 2.
K[i][k'][j][l] = The rate of flow of item 'i' from compartment 'k'', to compartment 'j'
\\Both Ks are time-independent
r[i][j][k'][l] = some flow coefficient same as K- which is also independent of time.
A[i][k'][j][l] = B[i][k'][j][l] = Process of the flow- can be either 0, 1 or 2.
I is also a constant and does not explicitly depend on time.

An example is
dx[i][j]/dt = I - K[i][j][k'][2]*pow(x[i][j],2)...

Which should answer the 'n' @lastchance asked. I should have clarified it isn't the same n as i had in my (t+nh).

2. https://imgur.com/a/4ksj7

This second image shows what I am trying to do with this model of mine. If we took model E for example, assuming there can be five unique items/substances/quantities present in each of (1,2,3) - which are compartments, we can say x[5][3] is a set of concentrations of such item/substances/quantities. This is what I am trying to get with my model, although I have 4 compartments, and 5 items which make it x[5][4] in my case as the differential equation above shows.

Last edited on

closed account (48T7M4Gy)

Last edited on

lastchance (6980)

@Kloppite,

Please could you supply:
- (ideally) an online web reference to your one of your supervisor's (or whoever's) journal papers with some testable data;

- for the above, if possible, or otherwise something workable and consistent:
--- the number of species [maximum i]
--- the number of locations [maximum j]
--- the set of initial concentrations Xij
--- the set of source terms Iij
--- the set of transfer coefficients (Kijk, Aijk, Gijk, Bijk and same with kj reversed)
Input data files will do, as long as you give us a rough idea of the format.

Last edited on

closed account (48T7M4Gy)

Last edited on

lastchance (6980)

#include <iostream>
#include <fstream>
#include <iomanip>
#include <valarray>
#include <vector>
using namespace std;

// Function prototypes
void step( double dx, double &x, valarray<double> &Y );         // one Runge-Kutta step
valarray<double> F( double x, valarray<double> Y );             // the derivative function (multi-variable)
void readNums();                                                // read the numbers of species and locations
void readSource();                                              // read all the source data (Iij)
void readTransfer();                                            // read all the transfer data (Kijk etc
void readConc( valarray<double> &Y );                           // read initial concentrations
void writeData( double t, valarray<double> Y, bool header );    // write data (optionally with header)
int ijIndex( int i, int j );                                    // convert from i,j to 1-d index
void indexij( int n, int &i, int &j );                          // convert from 1-d index to i,j

//*******
// Global variables
//*******
struct SOURCE
{
   int i, j;                // the species number (i) and location (j)
   double I;                // the source term Iij
};
istream & operator >> ( istream &strm, SOURCE &s )       // overload the input operator for a SOURCE
{ 
   strm >> s.i >> s.j >> s.I;
   return strm;
}
vector<SOURCE> sourcelist;                             


struct TRANSFER 
{
   int i, j, k;             // the species number (i), 'from' location (j) and 'to' location (k)
   double Kijk, Kikj;       // main rate constants     
   double Aijk, Aikj;       // rate-constant exponents 
   double Gijk, Gikj;       // the inertial constant   
   double Bijk, Bikj;       // inertial-constant exponents 
};
istream & operator >> ( istream &strm, TRANSFER &t )       // overload the input operator for a TRANSFER
{ 
   strm >> t.i >> t.j >> t.k >> t.Kijk >> t.Kikj >> t.Aijk >> t.Aikj >> t.Gijk >> t.Gikj >> t.Bijk >> t.Bikj;
   return strm;
}
vector<TRANSFER> transferlist;     


int NUMSPECIES;   
int NUMLOCATIONS; 


double dt = 0.1;            // timestep
int nstep = 100;            // number of timesteps
int nprint = 10;            // frequency of output
//*************
// End of global variables 
//*************

//========

int main()
{
   readNums();                                             // Get NUMSPECIES, NUMLOCATIONS
   valarray<double> Y(NUMSPECIES * NUMLOCATIONS);          // Set up 1-d array to hold dependent variables
   readSource();                                           // Read all source terms Iij into sourcelist
   readTransfer();                                         // Read all transfer parts of the equations into transferlist
   readConc( Y );                                          // read the initial concentrations

   double t = 0;
   writeData( 0, Y, true );                                // Output the starting conditions
   for ( int n = 1; n <= nstep; n++ )                      // Solve the differential equation using nstep intervals
   {
      step( dt, t, Y );                                    // Single step of Runge Kutta
      if ( n % nprint == 0 ) writeData( t, Y, false );     // Output every nprint steps
   }
}

//========

void step( double dx, double &x, valarray<double> &Y )     // Single Runge-Kutta step
{
   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;
}

//========

// YOUR MAIN DERIVATIVE (the RHS) IS SET HERE ***
valarray<double> F( double x, valarray<double> Y ) 
{
   valarray<double> f( Y.size() );               

   f = 0;                              // Start with 0

   for ( SOURCE s : sourcelist )       // Add the source terms
   {
      int i = s.i;
      int j = s.j;
      int nij = ijIndex( i, j );
      f[nij] += s.I;
   }

   for ( TRANSFER t : transferlist )   // Add the transfer terms 
   {
      int i = t.i;
      int j = t.j;
      int k = t.k;
      int nij = ijIndex( i, j );
      int nik = ijIndex( i, k );

      f[nij] += t.Kikj * pow( Y[nik], t.Aikj ) / ( 1.0 + t.Gikj * pow( Y[nik], t.Bikj ) );
      f[nij] -= t.Kijk * pow( Y[nij], t.Aijk ) / ( 1.0 + t.Gijk * pow( Y[nij], t.Bijk ) );
   }

   return f;
}

//========

void readNums() 
{
   ifstream infile( "numbers.dat" );
   infile >> NUMSPECIES >> NUMLOCATIONS;                   
   infile.close();
   cout << "Number of species: " << NUMSPECIES << "    Number of locations: " << NUMLOCATIONS << '\n';
}

//========

void readSource()
{
   SOURCE s;

   ifstream infile( "source.dat" );
   while ( infile >> s ) sourcelist.push_back( s );
   infile.close();
   cout << "Number of source items read: " << sourcelist.size() << '\n';
}

//========

void readTransfer()                
{
   TRANSFER t;
   ifstream infile( "transfer.dat" );
   while ( infile >> t ) transferlist.push_back( t );
   infile.close();
   cout << "Number of transfer items read: " << transferlist.size() << '\n';
}

//========

void readConc( valarray<double> &Y ) 
{
   int i, j;
   double X;

   Y = 0;    // Anything NOT set in the file ... is assumed to be 0

   ifstream infile( "concentration.dat" );
   while ( infile >> i >> j >> X ) Y[ ijIndex( i, j ) ] = X; 
   infile.close();
}

//========

int ijIndex( int i, int j )
{
   return NUMLOCATIONS * ( i - 1 ) + ( j - 1 );
}

//========

void indexij( int n, int &i, int &j )
{
   i = 1 + n % NUMLOCATIONS;
   j = 1 + n / NUMLOCATIONS;
}

//========

void writeData( double t, valarray<double> Y, bool header )
{
   #define SPH << setw( 10 ) <<
   #define SPD << setw( 12 ) << fixed << setprecision( 4 ) <<

   if ( header )
   {
      cout SPH "  " << " t";
      for ( int i = 1; i <= NUMSPECIES; i++ )
      {
         for ( int j = 1; j <= NUMLOCATIONS; j++ ) cout SPH 'C' << i << j;
      }
      cout << '\n';
   }

   cout SPD t;
   for ( int i = 1; i <= NUMSPECIES; i++ )
   {
      for ( int j = 1; j <= NUMLOCATIONS; j++ ) cout SPD Y[ ijIndex( i, j ) ];
   }
   cout << '\n';
}

Number of species: 2    Number of locations: 2
Number of source items read: 4
Number of transfer items read: 4
           t         C11         C12         C21         C22
      0.0000      1.0000      2.0000      0.1000      0.2000
      1.0000      1.1577      2.3423      0.1832      0.2068
      2.0000      1.3150      2.6850      0.2392      0.2408
      3.0000      1.4732      3.0268      0.2879      0.2821
      4.0000      1.6330      3.3670      0.3344      0.3256
      5.0000      1.7947      3.7053      0.3804      0.3696
      6.0000      1.9584      4.0416      0.4263      0.4137
      7.0000      2.1240      4.3760      0.4723      0.4577
      8.0000      2.2916      4.7084      0.5185      0.5015
      9.0000      2.4610      5.0390      0.5649      0.5451
     10.0000      2.6321      5.3679      0.6115      0.5885

The source files are:

numbers.dat (just NUMSPECIES and NUMLOCATIONS)

2   2

source.dat (any non-zero source term Iij; see overloaded operator >> for struct SOURCE)

1   1    0.2
1   2    0.3
2   1    0.05
2   2    0.04

transfer.dat (any non-zero equation coefficients; see overloaded operator >> for struct TRANSFER)

1  1  2    0.4   0.4   1.0   1.0   2.5   2.5   2   2
1  2  1    0.4   0.4   1.0   1.0   2.5   2.5   2   2
2  1  2    0.7   0.7   1.0   1.0   1.5   1.5   2   2
2  2  1    0.7   0.7   1.0   1.0   1.5   1.5   2   2

concentration.dat (any non-zero initial concentrations; see routine readConc() )

Last edited on

Kloppite (30)

Hi guys, thanks for the ideas so far. I haven't so much of a reference on this sort of model other than experimental data sets which are pretty much of little use. @Lastchance

Here is a format of one of my input file for the i = 1 that I have written:

1. https://imgur.com/a/qe4z7

First column = flag(which is supposed to serve as a marker for the rows when they have been minimized (later in the model))

Second column: item-number(i=1) here

Third Column: current location number (These files are written in terms of locations so there are four of them). This particular one is for location j = 1.

Fourth column: Location flowing to(from 0 - 4). I have included Zero here, however, we don't need values for concentrations at zero. k'= {0...4}

Fifth is the kinetic process number 'n' = (0 to 2).

Sixth and Seventh: Are A and B..These two take the values of the indexes in the four previous columns that is A[1][j][k'][n] = B[1][j][k'][n] = n.

Eighth and Ninth: Are K and Gamma..These two take the values of the indexes in the four columns from above so that the last two columns of each row is of the form: K[1][j][k'][n], Gamma[1][j][k'][n].

These are some test values though. I plan to start it out with values 1, for the K and Gamma.

Since i wrote these set of parameter files in terms of positions(j), here is one with j = 2.
https://imgur.com/a/gUiYJ

j =3&4
https://imgur.com/a/NbHgd

As i have mentioned this is just for the case of i = 1 that i have written and tested. These files have been updated for the entire i = 5 now in each compartment file.

Here is my list of items, which i think @kemort may be interested in

x[1][1],x[1][2],x[1][3],x[1][4] i =1, j = {1..4} (Measured in all four compartments)

x[2][1],x[2][2],x[2][3] i = 2, j = {1..3} (Measured in three compartments)

x[3][2] i = 3, j = {2} (Measured in one compartment = 2)

x[4][2] i = 4, j = {2} (Measured in one compartment = 2)

x[5][1],x[5][2] i = 5, j = {1-2}(Measured in two compartments = 1,2)

The initial conditions of all these concentrations i.e x[i][j](t_0) = 1.0.

I am trying to solve Rk method for each of these concentrations within each compartment that they are being measured.

I[i][j] is only specific to item 5, and then it is implicitly dependent on time - even though it doesn't change in time.

--- the number of species [maximum i] = 5
--- the number of locations [maximum j] = 4

Hope this clears it up a bit.

closed account (48T7M4Gy)

@lastchance & @kloppite

As you can see I decided to remove the earlier posts because after a bit of searching around I don't think they are relevant to this problem except for my comments on the large number of array dimensions.

If I'm right and this is about pharmokinetics compartment modelling as the textbook chapter 1 exercise indicates then these three references are useful to expand on it, and also to simplify and structure the number crunching:

1. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2811643/#Sec8
2. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3150575/
3. https://www.jstage.jst.go.jp/article/cpb1958/25/6/25_6_1312/_article

The matrix approach in the first one handles multi-compartment (connectivity) systems, and convolution/deconvolution is referred to in detail in the third.

I assume Gamma is the gamma function and if so Gamma[1] is incorrect. In C++ Gamma[1] is the 2nd element in an array called Gamma, whereas in the equation it is Gamma(1) ie the Gamma function with a parameter of 1.
( https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3150575/ )

I guess the bottom line for me is that I would:
- reduce/eliminate the number of dimensions dramatically. Subscripts are not generally dimensions.
- review whether 10 data files are needed - 2 at most - 1 for input & 1 for output
- use arrays (or better <vector>'s)
- review whether Runge-Kutta is needed at all, but separately, integrate with Simpson's rule or even with trapezoidal which is probably good enough?
- consider structs (or classes ) to encapsulate the data at some (later) stage
- get it to work for a single compartment -> then two -> then ...

If I get time I might try out a 1 and/or two compartment model and post here (permanently). The references have sufficient values for constants to get some realistic numbers to crunch.

Cheers :)

lastchance (6980)

Hi @kemort,

All I can say is that my last code post attempts to program up - in good faith - the set of equations given by the OP in his/her link
https://imgur.com/a/rpUou

I agree with you in so many points:-
- many dimensions are definitely not needed; internally my code puts all the dependent variables in a one-d (val)array; all the non-zero equation coefficients are held in 1-d lists (as vectors), not 4-dimensional arrays;
- definitely should include structs or classes to encapsulate data (at the moment I've simply created SOURCE and TRANSFER to represent content supplied from outside and transferred between compartments);
- reduce the input files; at the moment I haven't gone far enough, but I think the modelling is at a developmental stage and likely to evolve; I didn't want to risk fixing a format for a single input file which became unworkable later; it's also the reason my code holds on to what, to me, is an uncomfortable number of global variables.

The only point of disagreement is whether Runge-Kutta is necessary and whether you could just integrate. The equations given by the OP are non-linear in the dependent variables and appear (to me) impossible to rearrange for straight integration (except in the cases where a lot of the constants are zero). In the links that you give, the equations seem to be linear and could reasonably be solved by Laplace transforms/deconvolution. My experience of rate equations would indeed be with linear equations ... it's just not what the OP seems to be using. Believe me, if you've ever tried turbulence modelling in fluid mechanics you'd be staggered at how unnecessarily complex your colleagues want to make model equations, especially when they aren't the ones tasked with solving them.

(Just as a matter of interest, if you were trying to solve dY/dx = F(x,Y) and F() happened to be independent of Y then a Runge-Kutta solver automatically gives you a Simpson's-rule integrator for free: check the final weighting.)

Note to the OP: in the equations that you give there is no coupling between species at all; the same i index appears as the first index in every term. As given, you might as well solve for each species as a separate problem. Also, I can't see where "process number" comes into the problem.

Second note to the OP: for any rate equation a typical timescale is 1/K where K is the rate constant; for a set of processes the system needs to resolve time 1/Kmax. So, given the accuracy of Runge-Kutta, a suitable timestep would be about one tenth of that. To test for "convergence" (what I would probably call "time-step independence") just run with a time step half the size (and double the number of steps) and see if it gives the same result.

Last edited on

closed account (48T7M4Gy)

Hi @lastchance,
You always act in good faith so rest assured your hard work in coding wasn’t being questioned.

I haven’t had a chance to look closely at your program yet but it shows that OP can trim things down.

structs and classes also lend themselves to a network/graph model whereby any number of compartments connected in any number of configurations could be handled, in a model not unlike a complex linked list, directed graph.

I’ll revisit Runge Kutta.

Cheers :)

Kloppite (30)

Hi @lastchance and @Kemort

I really appreciate your suggestions and insight. A couple of those papers @Kemort suggested are actually good reads as an insight to what exactly I am trying to achieve, although mine is a bit more exhaustive.

In response to you, @Kemort, I am working on this with a mentor of mine who has worked on the model before although with Mathematica but insists on the effectiveness of the RK methods - with the dimensions that I have used. Probably for a measure of completeness, but I think as a reference, the RK method does work.

The Gamma term is also not so much a Gamma function as it is a flow constant of the equation. What it does is actually condenses the equation into a sort of Michaelis Menten equation.

@Lastchance

There is a coupling in the equation actually, but not in the image I posted. The coupling is observed only between i = 3 and i = 4 to give i = 2 in this sense

x[1][2] + x[3][2] \rigthleftharpoons x[4][2] \rigtharrow x[2][2] + x[3][2]

Which is an enzyme catalyzed reaction with a complex x[4][2].

Re: Convergence: Yeah, within my code I have written out a while loop with a convergence condition, within which I have step-size halved till I reach a tolerable limit with respect to the difference between successive values of x[1][j] of the divided step-sizes.

Something like pow(2, count) - 1: Steps for each time_step
where count increases within the While loop itself. However, this comes at a significant cost in terms of execution time.

Last edited on

Kloppite (30)

@Kemort @Lastchance
Hi, guys what do you think of this convergence function? I have written this extra code fragment to get the convergence of the concentration values in time after calculating the values at (t+h), stored in X, and values at (t) stored in Y so:

X = X(t+nh)
Y = X(t)

So I have set convergence at 100, to begin with, although I do have specific convergence test for each of the 20 concentrations X. The MaximumValue function iterates through all the 20 convergences to get the maximum and store it in variable 'conv' which is tested in the while loop.

My result shows a convergence of the one variable but not the other (i am still testing for 2 of the concentrations now). Any ideas?

Thank you.

void convX(double h_step, double t_i, valarray<double>&X, valarray<double>&Y, valarray<double>&X_conv){
valarray<double> conv_X(X.size());
double conv = 100.0;
conv_X = 100.0;
X_conv = Y;
int count = 1;
while((conv > 0.000001) && (count < 22)){

    h_step = h_step*0.5;
    int step_num_conv = static_cast<int>(pow(2,count)); // Tells the number of Iteration to get the final value for each count value
    for (int i_conv = 1; i_conv < step_num_conv+1; i_conv++){
       t_step(h_step, t_i, X_conv); // RK4 calculation using the Variable X_conv to store the values 
    }
    for (int i = 1; i <= DRUGNUM; i++){
        for (int j = 1; j <= COMPNUM; j++){
            if (X_conv[ijIndex(i,j) > -0.000001]){ // Convergence Limit 
            conv_X[ijIndex(i,j)] = fabs(X_conv[ijIndex(i,j)] - X[ijIndex(i,j)]); // Convergence test
            }
            else
            {
            conv_X[ijIndex(i,j)] = 100.00; //Setting the convergence to the initial value when the new Value X_conv is below an expected value
            }
        }
    }
    X = X_conv; // Setting the X_conv to X(t+h)
    conv = maximumValue(conv_X);// Testing for the maximum convergence value to know which concentration to further test.
    count++;
}

}

double maximumValue(valarray<double> &conv_X)
{
     int size_Conv = conv_X.size( );  // Establish size of array
     double MaxVal = conv_X[1];       // Start with max = first element

     for(int i = 0; i<size_Conv + 1; i++)
     {
          if(conv_X[i] > MaxVal)
                MaxVal = conv_X[i];
     }
     return MaxVal;                // Return highest value in array
}

Last edited on

closed account (48T7M4Gy)

@Kloppite
I picked on the easy one and you might like to consider a couple of things about style - a sometimes very subjective game. One primary goal is to make the code self documenting and easy to read, especially for later code maintenance and updating. Lines wider than 80 characters are ... bad, because they stuff up being displayed and/or printed.

It is also worth pointing out here that using STL containers like <valarray> (I'd use <vectors> but that's me) opens up a whole C++ technology and functionality in manipulating the data in the container.

See: http://en.cppreference.com/w/cpp/numeric/valarray, max is just a one liner.

If you have valarray<double> maybe call it Array, hence the typdef.

#include <iostream>
#include <valarray>

typedef std::valarray<double> Array;
double getMaximum(Array &); // functions are active DO things

int main()
{
    Array convergence_x {5,77,8,-9}; // big self-documenting names
    
    std::cout << getMaximum(convergence_x) << '\n';
    
    return 0;
}

//  PRE: Given an Array
// POST: Return maximum value
double getMaximum(Array &aArray)
{
    double maximum = aArray[0]; // 1st element is [0]
    
    for(size_t i = 0; i < aArray.size(); i++) // see STL begin(), end() etc
    {
        if(aArray[i] > maximum)
            maximum = aArray[i];
    }
    return maximum;
}

77
Program ended with exit code: 0

lastchance (6980)

Last edited on

closed account (48T7M4Gy)

Oops, I forgot:

//  PRE: Given an Array
// POST: Return maximum value
double getMaximum(Array &aArray)
{
    return aArray.max();
}

or just
std::cout << convergence_x.max() << '\n';

or

1
2

double what_is_maximum;
what_is_maximum = convergence.max;

closed account (48T7M4Gy)

@Kloppite,

I've looked at this function from the same point of view as the earlier one I commented on. I haven't delved into the mysteries of what the convergence is actually doing because I haven't spent any time on it. I gather it's some sort of reality check to make sure the RK answers are converging to a sensible result rather than spinning off into the wilderness.

I'm not a FORTRAN programmer so I'm biased a bit. I've got nothing against FORTRAN but FORTRAN code is not generally up to date on how to write clear programs. As you know, it was originally written to avoid hexcode but we're way past that stage.

(I hope I'm not offending the last few people still wearing white dust coats carrying around boxes of Hollerith cards.)

FWIW I think your function should return an Array if at all possible rather than passing values by reference. Either that or 'litter' your addresses and function with const keyword. The original RK4 stuff originally posted used passing parameters by reference but I'd avoid that by using the proverbial input->process->output model where the function is a machine to feed with input values and spits out the stuff you want.

I can't fathom your indexing which you called ijindex(i,j) and I called ndx etc to fit on a line. Something tells me that's not good.

Of course, it's your call, and I can accept it that you think yours is clearer than mine.

// PRE: ???
// POST: ???
Array getConvergence(double h, double t, Array &X, Array &Y) // one less Array?
{
    Array conv_X(X.size()); // ??? return this Array from the function
    
    double convergence = 100.0;
    conv_X = 100.0; // ???
    X_conv = Y; // ???
    
    int count = 0; // conventional wisdom is start at zero where possible
    
    double convergence_limit = -0.000001;
    
    while((convergence > -convergence_limit) && (count < 22))
    {
        h = h * 0.5;
        int max_no_iterations = static_cast<int>(pow(2,count));
        
        for (int iteration = 0; iteration < max_no_iterations; iteration++)
        {
            t_step(h, t, X_conv);
        }
        
        for (int i = 1; i <= DRUGNUM; i++)
        {
            for (int j = 1; j <= COMPNUM; j++)
            {
                // if (X_conv[ijIndex(i,j) > -0.000001])doesn't make sense to me
                
                if (X_conv[ndx(i,j)] > convergence_limit) // ???? maybe this
                {
                    conv_X[ndx(i,j)] = fabs(X_conv[ndx(i,j)] - X[ndx(i,j)]);
                }
                else
                {
                    conv_X[ndx(i,j)] = 100.00;
                }
            }
        }
        X = X_conv;
        
        convergence = maximumValue(conv_X);
        
        count++;
    }
    
    return conv_X; // ??? return
}

lastchance (6980)

You shouldn't change the timestep mid-calculation. If you want to change the global timestep then change the main loop (e.g. as below).

Note that the main advantage of using valarray is that it will do whole-array calculations in a single line of code (yes, @kemort, just like all modern FORTRAN arrays do by default!). e.g.
change = abs( Y - Ylast ).max();

This is as for my previous example. It isn't particularly well chosen as Runge-Kutta is already very accurate with a large step size.

int main()
{
   readNums();                                             // Get NUMSPECIES, NUMLOCATIONS
   valarray<double> Y(NUMSPECIES * NUMLOCATIONS);          // Set up 1-d array to hold dependent variables
   readSource();                                           // Read all source terms Iij into sourcelist
   readTransfer();                                         // Read all transfer parts of the equations into transferlist
   readConc( Y );                                          // read the initial concentrations

   double t = 0;
   writeData( 0, Y, true );                                // Output the starting conditions


   double dt = 0.1;                                           // Timestep
   int nstep = 100, nprint = 20;                              // Number of timesteps and output frequency
   const int MAXSTEP = 2000000;                               
   const double TOLERANCE = 0.00001;                          
   valarray<double> Yfirst = Y, Ylast = Y;                    // To store the start and end values
   double change = 1.0;                                       // To store the maximum absolute change from previous run

   while ( nstep <= MAXSTEP && change > TOLERANCE )
   {
      t = 0;   Y = Yfirst;                                    // Reset to initial conditions
      for ( int n = 1; n <= nstep; n++ )                      
      {
         step( dt, t, Y );
         if ( n % nprint == 0 ) writeData( t, Y, false );     // Output 
      }
      change = abs( Y - Ylast ).max();                        // Easy with valarray
      cout << "nstep = " << nstep << "      dt = " << dt << "      change = " << scientific << change << '\n' << string(40,'*') << '\n';
      nstep *= 2;   nprint *= 2;   dt /= 2.0;                 // Half the timestep and adjust steps accordingly
      Ylast = Y;
   }
}

Number of species: 2    Number of locations: 2
Number of source items read: 4
Number of transfer items read: 4
           t         C11         C12         C21         C22
      0.0000      1.0000      2.0000      0.1000      0.2000
      2.0000      1.3150      2.6850      0.2392      0.2408
      4.0000      1.6330      3.3670      0.3344      0.3256
      6.0000      1.9584      4.0416      0.4263      0.4137
      8.0000      2.2916      4.7084      0.5185      0.5015
     10.0000      2.6321      5.3679      0.6115      0.5885
nstep = 100      dt = 0.1000      change = 3.3679e+000
****************************************
      2.0000      1.3150      2.6850      0.2392      0.2408
      4.0000      1.6330      3.3670      0.3344      0.3256
      6.0000      1.9584      4.0416      0.4263      0.4137
      8.0000      2.2916      4.7084      0.5185      0.5015
     10.0000      2.6321      5.3679      0.6115      0.5885
nstep = 200      dt = 0.0500      change = 1.4439e-010
****************************************

As @kemort shows, writing valarray<double> every time begins to pall, so you can use aliases like
typedef valarray<double> Array;
or
using Array = valarray<double>;
However, I'm not going back to change all those usages now.

Last edited on

closed account (48T7M4Gy)

@lastchance
Regarding the step size, there's an interesting section in Numerical Recipes called
17.2 Adaptive Stepsize Control for Runge-Kutta that might be relevant here.

"Implementation of adaptive stepsize control requires that the stepping algorithm signal information about its performance, most important, an estimate of its truncation error." NUMERICAL RECIPES, The Art of Scientific Computing, Third Edition 2007

I don't speak with expertise on this but I casually encountered it only a couple of days ago.

Your defence of FORTRAN is noted, continue.

lastchance (6980)

Thanks @kemort.

There is a variant of Runge-Kutta (Runge-Kutta-Fehlberg) which is similarly 4th(or 5th?)-order accurate and explicit ... but also feeds back an error estimate.
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
This could certainly be used for adaptive step-size control. For problems that I have used it with it is so accurate that theoretical truncation errors are smaller than floating-point round-off - after which the law of diminishing returns sets in with a vengeance.

In this problem, however, since none of the rate constants change it is difficult to see why the timestep should change (except possibly to increase as the problem heads toward steady-state).

Oh, I can defend Fortran whenever; it is the language I use at work. If you read Chapman's book: "Fortran 95/2003 for engineers and scientists (3rd ed.)" you may be surprised at how Fortran has advanced: certainly to those who previously coded in Fortran 77 (or earlier). However, it had little to do with the choice of valarray here. That was simply because I wanted to write things like A - B or abs( A ) and be able to rely on them doing whole-array operations, without having to write some form of loop over elements. I write vectors down on paper with nothing more than a squiggle underneath: it is nice to be able to do similarly (without the squiggle!) in code.

Last edited on

Kloppite (30)

Hi, guys, thanks for the suggestions. My reasoning for using the convergence function is that I observed in some values that I do not get the convergence as I should.


         Time                X11                X12      
   -1.00000000    1.00000000    1.00000000
   -0.90000000    0.90000000    1.10000000
   -0.80000000    0.80000000    1.20000000
   -0.70000000    0.70000000    1.30000000
   -0.60000000    0.60000000    1.40000000
   -0.50000000    0.50000000    1.50000000
   -0.40000000    0.40000000    1.60000000
   -0.30000000    0.30000000    1.70000000
   -0.20000000    0.20000000    1.80000000
   -0.10000000    0.10000000    1.90000000
   -0.00000000    0.00000000    2.00416667
    0.10000000    0.00000000    2.00416667
    0.20000000    0.00000000    2.00416667
    0.30000000    0.00000000    2.00416667
    0.40000000    0.00000000    2.00416667
    0.50000000    0.00000000    2.00416667
    0.60000000    0.00000000    2.00416667
    0.70000000    0.00000000    2.00416667
    0.80000000    0.00000000    2.00416667
    0.90000000    0.00000000    2.00416667
    1.00000000    0.00000000    2.00416667

Here is a simple result from a test with k112 = 1, and k121 = 0.

The differential equations for the case above


dx11/dt = -k112*H(x11)
dx12/dt = k112*H(x11)

HERE H is a heavy side function that is zero when

 x11 <= 0.00

. This prevents a negative value of the differential equation.

The expected value of x12 from Time (-0.0000) should be

2.00000

or in the worst case a numerical noise of the same order as the convergence i.e. x12 =

2.000001

- Which is far different from the values above hence the inclusion of the convergence which doesn't appear to be working either.

Last edited on

Pages: 123

C++

Forum

4th order Runge-Kutta method of vectors