C++

Forum

 
4th order Runge-Kutta method of vectors

Pages: 123
lastchance (6980)
@kloppite,
To jump from 1.9000 to 2.00416667 over a timestep of 0.1 requires a gradient F bigger than 1. That is not possible for the equation you have given. Turn off any convergence test and look for the bug in either your code or your source files. I didn't code a step function in any of my samples and a power of zero doesn't correspond to a step function either. Repeat: this has nothing to do with "convergence". If you turn that off then you must be returning an F value bigger than 1 somewhere. Put some debugging lines in your derivative function to find out where.

Not the problem here, but if you are basing your Heaviside step function on a floating-point variable that is "exactly" 0 then you are risking floating-point arithmetic giving completely different results if it is actually -1e-20 or +1e-20. It really isn't a good test.
Last edited on
Kloppite (30)
Hey, lastchance.

Whats your idea of the best way to implement step functions using c++. I have to use a particular numeric step function for all my calculations, although I have done so using a bunch of if else statements, I believe it is not efficient in the long term since with convergence I may be having lots of calculation of certain steps. I have thought of using a Dirac function but will appreciate your input. Thanks.

The function is of the form

H(x) = {0 if X <= 0.0 and 1 otherwise.
lastchance (6980)
Hello, @Kloppite
I'm perfectly well aware of what a step function is. It's not a good test function here because floating-point imprecision is such that being a tiny bit less than 0, exactly equal to 0 or a tiny bit bigger than 0 makes a huge difference to something with a discontinuity of 1.0 in it at the origin.

However, that is NOT the problem here. Somewhere you are returning a gradient bigger than 1.0. That is a bug in your code, irrespective of the step function.

If you want a simple test function just use simple exponential decay; in your notation:
dx11/dt = -k112*x11

This gives a perfectly acceptable positive solution:
x11 = initial value * exp( -k112 * t )

and has the advantage of (a) being physically realistic (a step function isn't in the present context) and (b) smooth.

But I go back to the original premise. You have a gradient bigger than 1 in your equation for the other variable x12. Therefore, you have a bug in your program.

I have thought of using a Dirac function...

Please don't. Enough said.
Kloppite (30)
Thanks for the response. I have figured out the bug, just that the step function is an integral part of the code as it helps with continuous flow. Take for instance the equation you cited.
f1(t,x) = dx11/dt = -k112*x11

This will keep running even when x11 < 0, unless there is a condition such that returns zero when this happens hence the need for the step function.
lastchance (6980)
That exponential decay function is ALWAYS POSITIVE. You do NOT need a step function. The only way that your approximate numerical solution could go negative is if the timestep was too large. However, if you take a timestep less than 1/k that will not happen.

If you need to ensure that concentrations don't go negative then put a limiter on them at the end of a timestep. However, for the simple exponential decay this should not happen.
Last edited on
Kloppite (30)
Apologies @lastchance, that was just an illustration albeit a stupid one from my end. I meant a constant function, say,

f(t,x11) = dx11/dt = -k112
f(t,x12) = dx12/dt = k112
...............
...............
x11 += h/6(f(t,x11) + .... + f(t+h, x11+h*f(x11+...)))
x12 += h/6(f(t,x12) + .... + f(t+h, x12+h*f(x12+...)))


In this case, the limiter works for x11 but not x12, unless i add a H(x11){or H(x11+fx11)} function.
Last edited on
closed account (48T7M4Gy)
@Kloppite

Would you like to post your (10 or 16) equations in algebraic form with appropriate initial values and values for any constants? If you do, please include any restrictions on values eg the variables that must stay positive.
lastchance (6980)
@Kloppite, if you use an equation like
dx11/dt = -k112

then you can just write down the algebraic solution 
x11 = (initial value of x11) - k112.t


If you limit it via
dx11/dt = -k112 H(x11)

then, again, you can just write down the algebraic solution immediately:
x11 =max{  (initial value of x11) - k112.t ,  0 }


It's barely a differential equation. If you are going to use a complicated code to solve such as this then you are making a mountain out of a molehill.

I fail to see how "using a step function" could "help with continuous flow".
Last edited on
Kloppite (30)
@LastChance and Kemort

At Kemort's request, this is a form of my equations

I will use four differential equations here to represent 1 substances and 2 compartments.

I will assume the substances flow between compartments 1 and 2. So there is a presence of both drugs in the compartments, say, x11, x12.

The parameters still take the forms:
K{l}ijk = rate of flow of item 'i' in compartment 'j' to 'k' following process 'l'.
K{1}112 == rate of flow of item 1, in compartment 1 to compartment 2 using process 1.

The processes describe rate laws which can be Zeroth(0), First(1) or Second(2).

I is the source term, only flowing into compartment 1 and is a function of time, written I(t). In my current test it starts, it starts when (1.0>=t>=2.0).


Midpoint F1(f1,d1)

f1(t,x11) = I(t)-K{0}112*(x11)^{0} - K{1}112*(x11)^{1} - K{2}112*(x11)^{2} 
f1(t,x11) += K{0}121*(x12)^{0} + K{1}121*(x12)^{1} + K{2}121*(x12)^{2} 

d1(t,x12) =  -K{0}121*(x12)^{0} - K{1}121*(x12)^{1} - K{2}121*(x12)^{2}
d1(t,x12) +=K{0}112*(x11)^{0} + K{1}112*(x11)^{1} + K{2}112*(x11)^{2}


Midpoint F2(f2(t,x11+h*0.5*f1),d2(t,x12+h*0.5*d1))

f2(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f1)^{0} - K{1}112*(x11 + h*0.5*f1)^{1} - K{2}112*(x11 + h*0.5*f1)^{2} 

f2(t,x11) += K{0}121*(x12 + h*0.5*d1)^{0} + K{1}121*(x12 + h*0.5*d1)^{1} + K{2}121*(x12 + h*0.5*d1)^{2} 

d2(t,x12) = -K{0}121*(x12 + h*0.5*d1)^{0} - K{1}121*(x12 + h*0.5*d1)^{1} - K{2}121*(x12 + h*0.5*d1)^{2}

d2(t,x12) +=K{0}112*(x11 + h*0.5*f1)^{0} + K{1}112*(x11 + h*0.5*f1)^{1} + K{2}112*(x11 + h*0.5*f1)^{2}


Midpoint F3(f3(t,x11+h*0.5*f2),d3(t,x12+h*0.5*d2))

f3(t,x11) =  I(t+h/2) - K{0}112*(x11 + h*0.5*f2)^{0} - K{1}112*(x11 + h*0.5*f2)^{1} - K{2}112*(x11 + h*0.5*f2)^{2} 

f3(t,x11) += K{0}121*(x12 + h*0.5*d2)^{0} + K{1}121*(x12 + h*0.5*d2)^{1} + K{2}121*(x12 + h*0.5*d2)^{2} 

d3(t,x12) = -K{0}121*(x12 + h*0.5*d2)^{0} - K{1}121*(x12 + h*0.5*d2)^{1} - K{2}121*(x12 + h*0.5*d2)^{2}

d3(t,x12) +=K{0}112*(x11 + h*0.5*f2)^{0} + K{1}112*(x11 + h*0.5*f2)^{1} + K{2}112*(x11 + h*0.5*f2)^{2}


Midpoint F4(f4(t,x11+h*f3),d4(t,x12+h*d3))

f4(t,x11) = I(t+h) - K{0}112*(x11 + h*0.5*f3)^{0} - K{1}112*(x11 + h*0.5*f3)^{1} - K{2}112*(x11 + h*0.5*f3)^{2} 

f4(t,x11) += K{0}121*(x12 + h*d3)^{0} + K{1}121*(x12 + h*d3)^{1} + K{2}121*(x12 + h*d3)^{2} 

d4(t,x12) = -K{0}121*(x12 + h*d3)^{0} - K{1}121*(x12 + h*d3)^{1} - K{2}121*(x12 + h*d3)^{2}

d4(t,x12) +=K{0}112*(x11 + h*f3)^{0} + K{1}112*(x11 + h*f3)^{1} + K{2}112*(x11 + h*f3)^{2}


x11 += (h/6)*(f1 + 2f2 + 2f3 + f4)
x12 += (h/6)*(d1 + 2d2 + 2d3 + d4)

t starts at -1 to 5 minutes.





Last edited on
Kloppite (30)
In my tests, I am running specific cases and testing against theoretical values. My first case has 

I(t) = 15.0
K{0}112 = 1.0
x11 = x12 = 1.0

In addition to these, all other parameters are ZERO.

My equation now become
Midpoint F1(f1,d1)

f1(t,x11) = I(t)-K{0}112*(x11)^{0}

d1(t,x12) =K{0}112*(x11)^{0}


Midpoint F2(f2(t,x11+h*0.5*f1),d2(t,x12+h*0.5*d1))


f2(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f1)^{0}

d2(t,x12) =K{0}112*(x11 + h*0.5*f1)^{0}


Midpoint F3(f3(t,x11+h*0.5*f2),d3(t,x12+h*0.5*d2))


f3(t,x11) =  I(t+h/2) - K{0}112*(x11 + h*0.5*f2)^{0}

d3(t,x12) = K{0}112*(x11 + h*0.5*f2)^{0}


Midpoint F4(f4(t,x11+h*f3),d4(t,x12+h*d3))


f4(t,x11) = I(t+h) - K{0}112*(x11 + h*0.5*f3)^{0}

d4(t,x12) =K{0}112*(x11 + h*f3)^{0}




x11 += (h/6)*(f1 + 2f2 + 2f3 + f4)
x12 += (h/6)*(d1 + 2d2 + 2d3 + d4)


t runs from -1.0 to 5.0, and h = 0.1.

The concentrations are not allowed to be less than 10^{-6}, at which point they are set to ZERO.

xij > 10^{-6} ? xij : 0.0

My results from t = -0.9 to 1.0 (I == 0.0)

t                    X11                X12
-0.9          0.9000000      1.1000000
-0.80         0.8000000      1.2000000
-0.70         0.7000000      1.3000000
-0.60         0.6000000      1.4000000
-0.50         0.5000000      1.5000000
-0.40         0.4000000      1.6000000
-0.30         0.3000000      1.7000000
-0.20         0.2000000      1.8000000
-0.10         0.1000000      1.9000000
-0.00         0.0000000      2.0000000
0.10          0.0000000      2.0000000
0.20          0.0000000      2.0000000
0.30          0.0000000      2.0000000
0.40          0.0000000      2.0000000
0.50          0.0000000      2.0000000
0.60          0.0000000      2.0000000
0.70          0.0000000      2.0000000
0.80          0.0000000      2.0000000
0.90          0.0000000      2.0000000
1.00          0.0000000      2.0000000


My results from t = 1 to 2.0 (I == 15.0)

t                     X11               X12
1.10          1.4166667      2.0833333
1.20          2.8166667      2.1833333
1.30          4.2166667      2.2833333
1.40          5.6166667      2.3833333
1.50          7.0166667      2.4833333
1.60          8.4166667      2.5833333
1.70          9.8166667      2.6833333
1.80          11.2166667      2.7833333
1.90          12.6166667      2.8833333
2.00          14.0166667      2.9833333


Theoretical Values for same case t = -0.9 to 1.0 (I == 0.0)

-0.90   0.90000000    1.10000000
-0.80   0.80000000    1.20000000
-0.70   0.70000000    1.30000000
-0.60   0.60000000    1.40000000
-0.50   0.50000000    1.50000000
-0.40   0.40000000    1.60000000
-0.30   0.30000000    1.70000000
-0.20   0.20000000    1.80000000
-0.10   0.10000000    1.90000000
-0.00   0.00000000    2.00000000
0.10    0.00000000    2.00000000
0.20    0.00000000    2.00000000
0.30    0.00000000    2.00000000
0.40    0.00000000    2.00000000
0.50    0.00000000    2.00000000
0.60    0.00000000    2.00000000
0.70    0.00000000    2.00000000
0.80    0.00000000    2.00000000
0.90    0.00000000    2.00000000
1.00    0.00000000    2.00000000


Theoretical Values for same case t = 1.0 to 2.0 (I == 15.0)

1.10    1.40000000    2.10000000
1.20    2.80000000    2.20000000
1.30    4.20000000    2.30000000
1.40    5.60000000    2.40000000
1.50    7.00000000    2.50000000
1.60    8.40000000    2.60000000
1.70    9.80000000    2.70000000
1.80    11.20000000    2.80000000
1.90    12.60000000    2.90000000
2.00    14.00000000    3.00000000


I guess this explains my need for convergence. However i am open to suggestions if you see something i am not seeing. Thanks.
Last edited on
lastchance (6980)
@kloppite,
Please state the DIFFERENTIAL EQUATIONS that you wish solved. Please do NOT confuse the issue with the midpoint numerical method.

If you wish to solve FOUR simultaneous equations then write down those FOUR equations (not just two of them) PLAINLY in the form
d/dt(x11) = ...
d/dt(x12) = ...
d/dt(x21) = ...
d/dt(x22) = ...

State the numerical values of ALL the parameters Kij, Iij etc.

State the initial values of all of x11, x12, x21, x22 or your problem is not well posed.


Your problem is the step function. It has a massive discontinuity right at the end of a time interval. Since the Runge-Kutta method uses a 1/6 contribution from either end of an interval it doesn't really know whether to take 0 or 1 here: it even depends on whether you use >= 0 or > 0 in the definition of the step function. I REPEAT: THIS IS A TERRIBLE WAY TO TEST YOUR CODE: IT IS AN ISSUE WITH A DISCONTINUITY, NOT WITH CONVERGENCE.

The following is what I think the code will give you for your test problem (as best as I can interprete it) if you use a tiny timestep. The problem is that it requires a silly timestep because it is trying to round off the end of the interval.

Number of species: 1    Number of locations: 2
Number of source items read: 1
Number of transfer items read: 2
           t         C11         C12
****************************************
  --- Lines omitted ---
****************************************
           t         C11         C12
     -1.0000      1.0000      1.0000
     -0.9000      0.9000      1.1000
     -0.8000      0.8000      1.2000
     -0.7000      0.7000      1.3000
     -0.6000      0.6000      1.4000
     -0.5000      0.5000      1.5000
     -0.4000      0.4000      1.6000
     -0.3000      0.3000      1.7000
     -0.2000      0.2000      1.8000
     -0.1000      0.1000      1.9000
      0.0000      0.0000      2.0000
      0.1000      0.0000      2.0000
      0.2000      0.0000      2.0000
      0.3000      0.0000      2.0000
      0.4000      0.0000      2.0000
      0.5000      0.0000      2.0000
      0.6000      0.0000      2.0000
      0.7000      0.0000      2.0000
      0.8000      0.0000      2.0000
      0.9000      0.0000      2.0000
      1.0000      0.0001      2.0000
      1.1000      1.4001      2.1000
      1.2000      2.8001      2.2000
      1.3000      4.2001      2.3000
      1.4000      5.6001      2.4000
      1.5000      7.0001      2.5000
      1.6000      8.4001      2.6000
      1.7000      9.8001      2.7000
      1.8000     11.2001      2.8000
      1.9000     12.6001      2.9000
      2.0000     14.0001      3.0000
nstep: 120000   dt: 2.5000e-005   change: 2.8750e-004
****************************************

Last edited on
lastchance (6980)
Revised code - results to follow in a second post.

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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
 
#include <iostream>
#include <fstream>
#include <iomanip>
#include <valarray>
#include <vector>
#include <string>
using namespace std;

using Array = valarray<double>; 


//***** Function prototypes

void step( double dx, double &x, Array &Y );          // One Runge-Kutta step
valarray<double> F( double x, Array Y );              // Derivative function
void readNums();                                      // Read number of species and locations
void readSource();                                    // Read source data (Iij)
void readTransfer();                                  // Read transfer data (Kijk etc
void readConc( double &tfirst, Array &Y );            // Read initial concentrations
void writeData( double t, Array Y, bool header );     // Write data
int ijIndex( int i, int j );                          // Convert from i,j to n
void indexij( int n, int &i, int &j );                // Convert from n to i,j
double Apow( double x, int p );                       // Quasi-power law
void clip( Array &Y, double value );                  // Clip 


//***** Global variables

struct SOURCE
{
   int i, j;                // Species number (i) and location (j)
   double I;                // Source term Iij
   double tmin, tmax;       // Start & end times for Iij
};
istream & operator >> ( istream &strm, SOURCE &s )
{ 
   strm >> s.i >> s.j >> s.I >> s.tmin >> s.tmax;
   return strm;
}
vector<SOURCE> sourcelist;                             


struct TRANSFER 
{
   int i, j, k;             // Species number (i), 'from' location (j) and 'to' location (k)
   double Kijk, Kikj;       // Rate constants
   int Aijk, Aikj;          // Rate-constant exponents
   double Gijk, Gikj;       // Inertial constant
   int Bijk, Bikj;          // Inertial-constant exponents
};
istream & operator >> ( istream &strm, TRANSFER &t )
{ 
   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; 


//**** End of global variables


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

   double dt = 0.1;                                        // Timestep
   int nstep = 30, nprint = 1;                             // Initial number of timesteps and output frequency
   const int MAXSTEP = 200000;
   const double TOLERANCE = 0.00001;                       

   Array Yfirst = Y, Ylast = Y;                            // To store the start and end values
   double change = 1.0;                                    // To store the maximum absolute change from previous run
   double t;

   while ( nstep <= MAXSTEP && change > TOLERANCE )
   {
      t = tfirst;   Y = Yfirst;                            
      writeData( t, Y, true );                             // Output the starting conditions
      for ( int n = 1; n <= nstep; n++ )                   
      {
         step( dt, t, Y );                                 // One timestep
         clip( Y, 0.0 );                                   // Clip if necessary
         if ( n % nprint == 0 ) writeData( t, Y, false );  // Output 
      }
      change = abs( Y - Ylast ).max();                     
      cout << "nstep: " << nstep << scientific << "   dt: " << dt << "   change: " << change << '\n' << string(40,'*') << '\n';

      nstep *= 2;   nprint *= 2;   dt /= 2.0;              // Half the global timestep before repeat
      Ylast = Y;
   }
}

//========

void step( double dx, double &x, Array &Y )                // Single timestep
{
   int N = Y.size();
   Array dY1(N), dY2(N), dY3(N), dY4(N);

   // Runge-Kutta
   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;

   // Mid-point method
// dY1 = F( x           , Y             ) * dx;
// dY2 = F( x + 0.5 * dx, Y + 0.5 * dY1 ) * dx;
// Y += dY2;

   x += dx;
}

//========

Array F( double x, Array Y )                               // Derivative
{
   Array 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 );
      if ( x > s.tmin && x <= s.tmax ) 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 * Apow( Y[nik], t.Aikj ) / ( 1.0 + t.Gikj * Apow( Y[nik], t.Bikj ) );
      f[nij] -= t.Kijk * Apow( Y[nij], t.Aijk ) / ( 1.0 + t.Gijk * Apow( 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( double &tfirst, Array &Y )
{
   int i, j;
   double X;

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

   ifstream infile( "concentration.dat" );
   infile >> tfirst;
   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, Array 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';
}

//========

double Apow( double x, int p )                   // QUASI-power law
{
   switch ( p )
   {
      case 0: return ( x >= 1.0e-10 ? 1.0 : 0.0 );   // Step function
      case 1: return x;                              // Linear
      case 2: return x * x;                          // Quadratic
      default: 
         cout << "Faulty exponent in Apow";
         exit( 1 );
   }
}

//========

void clip( Array &Y, double value )
{
   int N = Y.size();
   for ( int i = 0; i < N; i++ ) if ( Y[i] < value ) Y[i] = value;
}

//========

lastchance (6980)
The code has been set up for the general set of equations you presented at
https://imgur.com/a/rpUou

So what does this code give for your extremely reduced "test problem", which I'm guessing is:
d(x11)/dt = I11(t) - K112 H(x11)
d(x12)/dt = K112 H(x11)     (<=== more about this monstrosity of an equation below)


Boundary conditions:
x11 = x12 = 1 at t = -1


Source: I11 = 15.0 if ( 1.0 < t <= 2.0 ) and 0.0 otherwise.
Rate constant: K112 = 1.0



The input files are:
numbers.dat: 
1   2


concentration.dat: - starts at time t = -1.0, with x11=1.0, x12=1.0
-1.0
1  1  1.0
1  2  1.0


transfer.dat: - massive overkill for this problem
1  1  2    1.0   0.0   0     0     0.0   0.0     0     0
1  2  1    0.0   1.0   0     0     0.0   0.0     0     0


source.dat: - source I11 = 15.0 from time 1.0 to 2.0 only
1   1    15.0    1.0   2.0



So what does this give with Runge-Kutta and some automatic global timestep convergence"?

Number of species: 1    Number of locations: 2
Number of source items read: 1
Number of transfer items read: 2
           t         C11         C12
     -1.0000      1.0000      1.0000
     -0.9000      0.9000      1.1000
     -0.8000      0.8000      1.2000
     -0.7000      0.7000      1.3000
     -0.6000      0.6000      1.4000
     -0.5000      0.5000      1.5000
     -0.4000      0.4000      1.6000
     -0.3000      0.3000      1.7000
     -0.2000      0.2000      1.8000
     -0.1000      0.1000      1.9000
     -0.0000      0.0167      1.9833
      0.1000      0.0000      2.0333
      0.2000      0.0000      2.0333
      0.3000      0.0000      2.0333
      0.4000      0.0000      2.0333
      0.5000      0.0000      2.0333
      0.6000      0.0000      2.0333
      0.7000      0.0000      2.0333
      0.8000      0.0000      2.0333
      0.9000      0.0000      2.0333
      1.0000      0.0000      2.0333
      1.1000      1.2000      2.0833
      1.2000      2.6000      2.1833
      1.3000      4.0000      2.2833
      1.4000      5.4000      2.3833
      1.5000      6.8000      2.4833
      1.6000      8.2000      2.5833
      1.7000      9.6000      2.6833
      1.8000     11.0000      2.7833
      1.9000     12.4000      2.8833
      2.0000     13.5500      2.9833
nstep: 30   dt: 1.0000e-001   change: 1.2550e+001
****************************************

 ----- LOTS OF LINES OMITTED -----

****************************************
           t         C11         C12
     -1.0000      1.0000      1.0000
     -0.9000      0.9000      1.1000
     -0.8000      0.8000      1.2000
     -0.7000      0.7000      1.3000
     -0.6000      0.6000      1.4000
     -0.5000      0.5000      1.5000
     -0.4000      0.4000      1.6000
     -0.3000      0.3000      1.7000
     -0.2000      0.2000      1.8000
     -0.1000      0.1000      1.9000
     -0.0000      0.0000      2.0000
      0.1000      0.0000      2.0000
      0.2000      0.0000      2.0000
      0.3000      0.0000      2.0000
      0.4000      0.0000      2.0000
      0.5000      0.0000      2.0000
      0.6000      0.0000      2.0000
      0.7000      0.0000      2.0000
      0.8000      0.0000      2.0000
      0.9000      0.0000      2.0000
      1.0000      0.0000      2.0000
      1.1000      1.4000      2.1000
      1.2000      2.8000      2.2000
      1.3000      4.2000      2.3000
      1.4000      5.6000      2.4000
      1.5000      7.0000      2.5000
      1.6000      8.4000      2.6000
      1.7000      9.8000      2.7000
      1.8000     11.2000      2.8000
      1.9000     12.6000      2.9000
      2.0000     13.9999      3.0000
nstep: 122880   dt: 2.4414e-005   change: 1.2207e-005
****************************************



The number of steps and the eventual timestep are ridiculous for this problem.

The reason is that you have a DISCONTINUITY in the derivative right at the point where Runge-Kutta is trying to get a value. With a step function H(x) there is NO DEFINED VALUE at x = 0: you can define (or fudge) it in any way you choose. Certainly H(x) = 0 for x < 0, and H(x) = 1 for x > 0 ... but what value should it have when x is precisely 0?

Then look at your second equation:
d(x12)/dt = K112 H(x11)


Depending on whether you define H(0) = 0 or 1, when x11 reaches 0 then this is EITHER
d(x12)/dt = 0
OR
d(x12)/dt = 1

These are completely different equations: that is why your test case is ambiguous and silly.



As a further consideration look inside the step() function: comment out the RUNGE-KUTTA lines and un-comment the MID-POINT method to use that instead. This time your output is
Number of species: 1    Number of locations: 2
Number of source items read: 1
Number of transfer items read: 2

--- LINES OMITTED ---

****************************************
           t         C11         C12
     -1.0000      1.0000      1.0000
     -0.9000      0.9000      1.1000
     -0.8000      0.8000      1.2000
     -0.7000      0.7000      1.3000
     -0.6000      0.6000      1.4000
     -0.5000      0.5000      1.5000
     -0.4000      0.4000      1.6000
     -0.3000      0.3000      1.7000
     -0.2000      0.2000      1.8000
     -0.1000      0.1000      1.9000
      0.0000      0.0000      2.0000
      0.1000      0.0000      2.0000
      0.2000      0.0000      2.0000
      0.3000      0.0000      2.0000
      0.4000      0.0000      2.0000
      0.5000      0.0000      2.0000
      0.6000      0.0000      2.0000
      0.7000      0.0000      2.0000
      0.8000      0.0000      2.0000
      0.9000      0.0000      2.0000
      1.0000      0.0000      2.0000
      1.1000      1.4000      2.1000
      1.2000      2.8000      2.2000
      1.3000      4.2000      2.3000
      1.4000      5.6000      2.4000
      1.5000      7.0000      2.5000
      1.6000      8.4000      2.6000
      1.7000      9.8000      2.7000
      1.8000     11.2000      2.8000
      1.9000     12.6000      2.9000
      2.0000     14.0000      3.0000
nstep: 120   dt: 2.5000e-002   change: 4.4409e-015
****************************************


This time the solution is obtained with only a small number of steps and a sensible timestep.

But this does NOT mean that the mid-point rule is better than Runge-Kutta ... it is just that the midpoint rule (starting at t = -1.0) DOES NOT USE A VALUE RIGHT ON THE DISCONTINUITY. If you started at a different time point then you will hit exactly the same problems as before.

There are similar problems accruing from the discontinuity in the source term.

The only way that you are going to handle a discontinuity correctly is to stop your code at the point of discontinuity, reset the initial conditions for variables and change the values in the input files so that you have no step functions in your derivative function. Then restart. You could, in principle, do this in code also, but it would be messy and error-prone, extremely ad-hoc, very problem-dependent and completely unphysical.

Use a better test case.

To see just how discontinuous this test problem is ...
https://imgur.com/RoTXfye
Last edited on
closed account (48T7M4Gy)
@Kloppite, thanks for the info but I can't fathom most of it.

From a a non-expert and FWIW:
The two compartment pharmacokinetic system can be described by the two differential equations below - I understand the equations I have used are for a bolus injection and I know that an additional ODE can accommodate more complex drug delivery.

There are two components no. 1 and no. 2 with corresponding drug amounts A1 and A2, and flow constants between the components being Kxy, where xy indicates the (connectivity direction between compartment x and compartment y).

The conclusions I have reached and if I'm right are that all the subscripts @Kloppite has used other than the compartment number are a confusing waste of time. Array dimensions aren't needed. The equations are pretty simple, even if there were ten compartments. It's just not worth the trouble unless you wanted to automate the process to get machine written ODE's - then connectivity matrices would be a part of the solution.

The program below is just using euler integration and if you plot the results using a spreadsheet the flow between the two components follows two nice exponential curves to a steady state - A2 rises as the A1 curve falls. RK-4 probably isn't needed because by the look of it all it will do is smooth an already smooth curve. Keep in mind I just picked arbitrary constants. They were first pass and haven't been adjusted - the curves are "very powerful" and "tolerant" at decay)

There is no need with this to consider -ve values. They occur as @lastchance said as an artifact of the numbers/time slice value. If they go negative it's an indicator of a blooper I'd say.

I'm making a whole lot of assumptions, particularly whether this is the same problem @Kloppite is trying to solve. Use or discard as you see fit :)

It's interesting though, they're nice plots, and a few more lines would show elimination etc etc. The big challenge is to calibrate it against reality - time as one obvious dimension.

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
 
#include <iostream>
#include <iomanip>

#include <vector>
#include <functional>

typedef std::function<double(double, double, double)> ODE;
typedef std::vector<ODE> ODE_vector;

// ODE's
double dA1(double, double, double);
double dA2(double, double, double);

// FLOW RATE CONSTANTS
const double k10 = 3.0;
const double k12 = 2.0;
const double k21 = 2.5;

// TIME INCREMENT
const double h = 0.05;

int main()
{
    ODE_vector fn;
    fn.push_back(dA1);
    fn.push_back(dA2);
    
    double t = 0, t0 = 0;
    double y = 0, y0 = 0;
    
    t = t0;
    
    double A1 = 0.1;
    double A2 = 0.0;
    
    double A1_next = A1;
    double A2_next = A2;
    
    std::cout << std::fixed;
    
    std::cout
    << "    t         A1      A2\n"
    << "------------------------\n";
    
    for(int i = 0; i <= 20; i++)
    {
        std::cout
        << std::setprecision(2)
        << std::setw(5) << t << '\t'
        
        << std::setprecision(4)
        << std::setw(8) << A1_next
        << std::setw(8) << A2_next
        << '\n';
        
        A1_next = A1_next + h*fn[0](t, A1_next, A2_next); // EULER
        A2_next = A2_next + h*fn[1](t, A1_next, A2_next); // EULER
        
        t += h;
    }
}

// COMPARTMENT 1 - CENTRAL
double dA1(double t, double A1, double A2)
{
    return k21*A2 - k12*A1 - k10*A1;
}

// COMPARTMENT 2 - PERIPHERAL
double dA2(double t, double A1, double A2)
{
    return k12*A1 - k21*A2;
}


 t         A1      A2
------------------------
 0.00	  0.1000  0.0000
 0.05	  0.0750  0.0075
 0.10	  0.0572  0.0123
 0.15	  0.0444  0.0152
 0.20	  0.0352  0.0168
 0.25	  0.0285  0.0176
 0.30	  0.0236  0.0177
 0.35	  0.0199  0.0175
 0.40	  0.0171  0.0170
 0.45	  0.0150  0.0164
 0.50	  0.0133  0.0157
 0.55	  0.0119  0.0149
 0.60	  0.0108  0.0141
 0.65	  0.0099  0.0133
 0.70	  0.0091  0.0126
 0.75	  0.0084  0.0118
 0.80	  0.0078  0.0111
 0.85	  0.0072  0.0105
 0.90	  0.0067  0.0098
 0.95	  0.0063  0.0092
 1.00	  0.0059  0.0087
Program ended with exit code: 0

https://imgur.com/gallery/0qRno
Last edited on
lastchance (6980)
particularly whether this is the same problem @Kloppite is trying to solve


@Kloppite's test problem solution is plotted at
https://imgur.com/RoTXfye
Last edited on
closed account (48T7M4Gy)
That display could be a dead patient or a Nobel prize - I'll have to wait for @Kloppite. The time-warp of -ve time is interesting - could be a combined medicine-physics-chemistry gong coming up :)
lastchance (6980)
The time-warp of -ve time is interesting


Sadly, that's what he wanted - mine is not to reason why.
Topic archived. No new replies allowed.
Pages: 123