0

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