Runge Kutta with C

[QuantuMechanic]

Greetings all,

I have written the following code to calculate the solution to a system of ODEs, called the Matsuoka equations, by using the Runge-Kutta 4th order method. I am trying to implement this code on an Arduino microcontroller. The output of the equations, IC[0] - IC[2], should oscillate but instead it "blows up (down)" and goes to negative infinity. This is written for Arduino so some things may need to be changed to get it to work in a pure c or c++ environment.

First I set up 5 arrays: 1 dummy array and 4 arrays to hold the values of the K1,K2,K3, and K4 coefficients. The Runge Kutta algorithm is set up as a function. The program calls the function and passes the dummy K array by reference to the function. The K array is altered due to the function calculations and then those values are stored in the K1 array to be used later to calculate the next step of the ODE. This process is repeated 3 times for the K2,K3, and K4 arrays. Finally the next step of the solution to the ODE is solved using K1,K2,K3, and K4.

My problem is made curiouser by the fact that I have gotten the code to work by foregoing all the for loops and solving each K coefficient line-by-line. This code is displayed below the NON-WORKING CODE and is termed WORKING CODE. What's interesting here is that the two codes produce the same values for the K1 and K2 arrays but for the K3 array the values differ in the K3[4] array position. The working code calculates K3[4] = 0.5 (which I confirmed by calculation) while the non-working code calculates K3[4] = 0.49999995230. Why does it round?

Any help as to why the NON-WORKING CODE blows up while the WORKING CODE oscillates would be greatly apperciated.

NON-WORKING CODE

void setup () {
Serial.begin(9600);
}

float c = 1, B = 2.5, Y = 2.5, tau1 = 1, tau2 = 2*tau1, timestep = 0.1, M = 0.5*timestep;
float IC[] = {1,1,1,0}, K[] = {0,0,0,0}, K1[3], K2[3], K3[3], K4[3];
int i;

void loop() {

rungeKutta(tau1,tau2,c,B,Y,0,IC[0],IC[1],IC[2],IC[3],K); //note M has been replaced by zero
for(i=0;i<=3;i++)
K1 = K;

rungeKutta(tau1,tau2,c,B,Y,M,IC[0],IC[1],IC[2],IC[3],K);
for(i=0;i<=3;i++)
K2 = K;

rungeKutta(tau1,tau2,c,B,Y,M,IC[0],IC[1],IC[2],IC[3],K);
for(i=0;i<=3;i++)
K3 = K;

rungeKutta(tau1,tau2,c,B,Y,timestep,IC[0],IC[1],IC[2],IC[3],K); //note M has been replaced by timestep
for(i=0;i<=3;i++)
K4 = K;

for(i=0;i<=3;i++)
IC = IC + (1.0/6.0)*(K1 + 2.0*(K2 + K3) + K4)*timestep;

Serial.println(IC[0] - IC[2],DEC);

}//end main loop

float rungeKutta(float tau1, float tau2, float c, float B, float Y, float M, float IC0, float IC1, float IC2, float IC3, double *K)
{
K[0] = (1.0/tau1)*(c - (IC0 + M*K[0]) - B*(IC1 + M*K[0]) - Y*max(IC2 + M*K[0],0));
K[1] = (1.0/tau2)*(max(IC0 + M*K[1],0) - IC1 + M*K[1]);
K[2] = (1.0/tau1)*(c - (IC2 + M*K[2]) - (B*IC3 + M*K[2]) - Y*max(IC0 + M*K[2],0));
K[3] = (1.0/tau2)*(max(IC2 + M*K[3],0) - (IC3 + M*K[3]));
}//end of rungeKutta function

WORKING CODE

void setup () {

Serial.begin(9600);

}

float c = 1, B = 2.5, Y = 2.5, tau1 = 1, tau2 = 2*tau1, timestep = 0.1, M = 0.5*timestep;;
float IC[] = {1,1,1,0}, K1[3], K2[3], K3[3], K4[3];

void loop() {

K1[0] = (1.0/tau1)*(c - IC[0] - B*IC[1] - Y*max(IC[2],0));
K1[1] = (1.0/tau2)*(max(IC[0],0) - IC[1]);
K1[2] = (1.0/tau1)*(c - IC[2] - B*IC[3] - Y*max(IC[0],0));
K1[3] = (1.0/tau2)*(max(IC[2],0) - IC[3]);

K2[0] = (1.0/tau1)*(c - (IC[0] + M*K1[0]) - B*(IC[1] + M*K1[0]) - Y*max(IC[2] + M*K1[0],0));
K2[1] = (1.0/tau2)*(max(IC[0] + M*K1[1],0) - IC[1] + M*K1[1]);
K2[2] = (1.0/tau1)*(c - (IC[2] + M*K1[2]) - (B*IC[3] + M*K1[2]) - Y*max(IC[0] + M*K1[2],0));
K2[3] = (1.0/tau2)*(max(IC[2] + M*K1[3],0) - (IC[3] + M*K1[3]));

K3[0] = (1.0/tau1)*(c - (IC[0] + M*K2[0]) - B*(IC[1] + M*K2[0]) - Y*max(IC[2] + M*K2[0],0));
K3[1] = (1.0/tau2)*(max(IC[0] + M*K2[1],0) - IC[1] + M*K2[1]);
K3[2] = (1.0/tau1)*(c - (IC[2] + M*K2[2]) - (B*IC[3] + M*K2[2]) - Y*max(IC[0] + M*K2[2],0));
K3[3] = (1.0/tau2)*(max(IC[2] + M*K2[3],0) - (IC[3] + M*K2[3]));

K4[0] = (1.0/tau1)*(c - (IC[0] + timestep*K3[0]) - B*(IC[1] + timestep*K3[0]) - Y*max(IC[2] + timestep*K3[0],0));
K4[1] = (1.0/tau2)*(max(IC[0] + timestep*K3[1],0) - IC[1] + timestep*K3[1]);
K4[2] = (1.0/tau1)*(c - (IC[2] + timestep*K3[2]) - (B*IC[3] + timestep*K3[2]) - Y*max(IC[0] + timestep*K3[2],0));
K4[3] = (1.0/tau2)*(max(IC[2] + timestep*K3[3],0) - (IC[3] + timestep*K3[3]));

IC[0] = IC[0] + (1.0/6.0)*(K1[0] + 2.0*K2[0] + 2.0*K3[0] + K4[0])*timestep;
IC[1] = IC[1] + (1.0/6.0)*(K1[1] + 2.0*K2[1] + 2.0*K3[1] + K4[1])*timestep;
IC[2] = IC[2] + (1.0/6.0)*(K1[2] + 2.0*K2[2] + 2.0*K3[2] + K4[2])*timestep;
IC[3] = IC[3] + (1.0/6.0)*(K1[3] + 2.0*K2[3] + 2.0*K3[3] + K4[3])*timestep;

Serial.println(IC[0] - IC[2],DEC);

}//end main loop

[mike_2000_17]

You need to make your arrays of size 4, not of size three. This means, your array declarations should be:

float IC[] = {1,1,1,0}, K[] = {0,0,0,0}, K1[4], K2[4], K3[4], K4[4]; //notice 4 here.

You basically has a problem of overlapping memory (the last element of K1 was at the same memory place as the first element of K2, as an example, but it could involve different elements). I'm pretty sure that would fix your problem.

>>K3[4] = 0.49999995230. Why does it round?
That's completely normal. You are using a "float" variable which is not expected to have better precision than to the 7-8th significant digit. This number is exactly what you would expect to get. You could use "double" for more precision (15 significant digits), but it's not sure that an Arduino uC would support that (sometimes floating point operations are not even supported at all).

[QuantuMechanic]

Thanks a ton Mike 2000, you are exactly right.