strange result

[genie0582]

Hi,

I use Runge Kutta to solve for ODEs. (Thanks for someone who provided this runge kutta code. I'm sorry if I didn't put a credit here)

I was wondering why when I set TIME 75, the result of y[3] becomes negative at the time around 70+. it is supposed to be zero as the same as when I use Mathematica.

Is there any modification on the equations before we use it? please help.

#include <stdio.h>
#include <stdlib.h>
#define TIME 75  
#define DELTA_TIME 0.1
#define N 4
#define dist DELTA_TIME
#define FI 20

double fp(double x, double y[],int k);

//Runge Kutta 4th order
void runge4(double x, double y[], double step)
{
    double H=(step/FI);
    double h=H/2.0,                      
        t1[N], t2[N], t3[N],            
        k1[N], k2[N], k3[N],k4[N];      
    int i,j;

for (j=1; j*H<=step ;j++)                     /* time loop */
{
    for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=H*fp(x, y, i));
    for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=H*fp(x+h, t1, i));
    for (i=0;i<N;i++) t3[i]=y[i]+    (k3[i]=H*fp(x+h, t2, i));
    for (i=0;i<N;i++) k4[i]=                H*fp(x+H, t3, i);
    for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
    x=j*H;
    }
}

int main(int argc, char *argv[])
{
    double t, y[N];
    int j;

    y[0]=50000.0; 
    y[1]=1.0;
    y[2]=1.0;
    y[3]=0.0;   
  
      for (j=1; j*dist<=TIME ;j++)
    {     
        runge4(t, y, dist);  
        t=j*dist; 
        printf("%.3lf   %.4lf   %.4lf   %.4lf    %.4lf \n",t,y[0],y[1],y[2],y[3]);
    } 
  
  system("PAUSE");	
  return 0;
}

//Kinetics
double fp(double x, double y[],int k) {
  switch(k) { 
    case 0: return -y[0];
    case 1: return 0.08144*(y[3]/(5000+y[3]))*y[1] - 0.005*y[1];
    case 2: return 0.1359*(y[3]/(0.05+y[3]))*y[2] - 0.009*y[2];
    case 3: return y[0] - 2.036*(y[3]/(5000+y[3]))*y[1] - 0.906*(y[3]/(0.05+y[3]))*y[2]; 
    default: return 0;
    }
}

[tesuji]

Hi,

I think your RK steps are not correct, you may check this.

Also a RK function of my own shows other steps for y direction:

tSCALAR ruku5( tODEF f, tSCALAR x, tSCALAR y, tSCALAR h)
// Common 4th order Runge–Kutta method
{  tSCALAR k1, k2, k3, k4;
   k1 = h * f ( x, y );
   k2 = h * f ( x + h / 2, y + k1 / 2 );
   k3 = h * f ( x + h / 2, y + k2 / 2 );
   k4 = h * f ( x + h    , y + k3 );
   return (y + ( k1 + 2 * ( k2 + k3 ) + k4 ) / 6);
}

So I think some of your loops should be changed to as commented:

for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=H*fp(x, y, i));
    for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=H*fp(x+h, t1, i));      // <-- H*fp(x+h,k1[i]/2, i)
    for (i=0;i<N;i++) t3[i]=y[i]+    (k3[i]=H*fp(x+h, t2, i));      // <-- H*fp(x+h, k2[i]/2,i)
    for (i=0;i<N;i++) k4[i]=                H*fp(x+H, t3, i);       // <-- H*fp(x+H, k3[i], i)
    for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;

Obviously the arrays t1, t2, t3 aren't necessary for computing y in last for-loop. Therefore scalar k1, k2, k3, k4, and array y can be computed within one loop.

-- tesu

[genie0582]

Thanks For your reply. I did what you suggested however I encountered another problem

for (j=1; j*H<=step ;j++)                     
{
     for (i=0;i<N;i++) k1[i]=H*fp(x, y, i);
     for (i=0;i<N;i++) k2[i]=H*fp(x+h, y[i]+0.5*k1[i], i);
     for (i=0;i<N;i++) k3[i]=H*fp(x+h, y[i]+0.5*k2[i], i);
     for (i=0;i<N;i++) k4[i]=H*fp(x+H, y[i]+    k3[i], i);

     for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
     x=j*H;
}

it reports "incompatible type for argument 2 of 'fp'
it means that we cannot just assign y+0.5*k1 in a function of fp(x,y,i).
That's why we still need t1, t2 and t3.

Any help?

[tesuji]

Sorry I have to apologize to you. At a first glance I thought your steps in y directions were wrong. But I was wrong. Your RK formulas are completely correct, as far as I can verify this, because of your system you need precomputing of the Ks and store them in t1..t3. Therfore the arrays are necessary. Can you post the results of Mathematica for time 70.3 and 75, like the result (t,y[0],y[1],y[2],y[3]) of printf?

Did you also consider that the results could be probably instable because 4th order RK (O(h^4)) might be insufficient?
-- tesu

[genie0582]

Thanks for your reply.
Your last assumption is what I was thinking also.
That's why I made FINESSE to the step size again.
however, it still gives negative values as shown in the table 1.
I was thinking maybe y[3] has a sharp turn at 71.4 - 71.5 (table 2)

If Runge Kutta 4th is not suitable for this problem. So what is the best solution? do you have any suggestion?
I'm thinking to use Runge Kutta 6th, but I believe the result will be the same.

Here are the results for both C language and Mathematica (NDSolve function)
Table 1
[IMG]http://farm5.static.flickr.com/4096/4742620854_08c88b6696.jpg[/IMG]
Sorry. I still don't know how to print the file from C.

Table 2
[IMG]http://farm5.static.flickr.com/4094/4742620924_2b2d6f6c78.jpg[/IMG]

[tesuji]

Hi

I have had a little time this morning so I did some tests with your c program. I run it using long double (12bytes with gnuc++) but the results were nearly the same as for double what is 8 byts long. Also I played about with some various step widths and got same bad results too.

Provided that the math, especially the kinetics, of your c program is exactly the same as for Mathematica, what you have most certainly verified, I almost feel certain that the procedural error of RK4 could be the reason for these bad results. Do you know which algorithm is implemented in Mathematica for solving ODE systems? (I myself are using matlab).

What do you think of Runge Kutta Fehlberg RK45 or RK4 with adaptive step size? Maybe adaptive step size leads to better results in time range 71.2 to 71.5 where the Yi feature strong alternations?

I wish you the best of luck

-- tesu

[genie0582]

Hi tesu,

I appreciate your kindness by checking on my code. Thanks again.
From the manual, built-in function in Mathematica (NDSOLVE) uses Method of Lines which are usually used to solve for PDEs in 1D. I don't have any idea about Method of Lines, so i cannot tell you much.

I have heard about RK45 with adaptive step size. When I check on the method, it looks so complicated. Anyway, I will try and let you know the result.

Thanks for you advice.

[genie0582]

I can figured it out now.

The ODEs I tried to solve are derived from biological model. This set of ODEs cannot be solved using any explicit scheme.
I posted for the help and then I helped myself.

[tesuji]

hi

glad to hearing of you. I was also afraid that the requirements (mainly initial value?) for integrating your partial diff. equations by means of a system of ODEs using RK solver weren't achieved. Unfortunately, I have lost so much of my knowledge on numerical how-to solve ODE/DAE/PDE ever since I left the uni for over ten years, so I wasn't any more able to give fair help.

I wish you all the very best for your research.

-- tesu