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
