I've been converting a Euler's method implementation of 4 linked differential equations, into a 4th order Runge Kutta implementation. I'm reasonably sure I've got the general approach right, and I've understood how to apply RK4, but I've not done any halfway serious maths for maybe 6 years so I may be missing something. My RK4 calculation gives sensible output when I use stepsize 1, but collapses to zero very quickly if I use any other step size. I'm hoping a fresh set of eyes might be able to quickly spot what I've done wrong. Please don't post a complete solution - I'd much prefer pointers to errors I might have made - whether that's code or understanding of RK4, as I'd like to be able to understand this for myself.
So here's my Euler implementation. Works well
// these defs are used in both Euler and RK4 versions
#define POSPART(X) (X > 0.0 ? X : 0.0)
#define NEGPART(X) (X < 0.0 ? X : 0.0)
#define STEP_DIVISION 10.0
double matsu_calc_nextVal(double in, double t1, double t2,
double c, double b, double g,
matsu_state *state)
{
double step = 1.0 / STEP_DIVISION;
double posX1, posX2, posIn, negIn, dx1, dx2, dv1, dv2;
double t1recip = 1.0/ t1;
double t2recip = 1.0/ t2;
for(int i=0; i<STEP_DIVISION; i++){
posX1 = POSPART(state->x1);
posX2 = POSPART(state->x2);
posIn = POSPART(in);
negIn = NEGPART(in);
dx1 = (c - state->x1 - (b*(state->v1)) - (posX2*g) - posIn) * t1recip;
dx2 = (c - state->x2 - (b*(state->v2)) - (posX1*g) - negIn) * t1recip;
dv1 = (posX1 - state->v1) * t2recip;
dv2 = (posX2 - state->v2) * t2recip;
state->x1 += dx1*step; state->x2 += dx2*step;
state->v1 += dv1*step; state->v2 += dv2*step;
}
return POSPART(state->x1) - POSPART(state->x2);
}
This is my RK4 implementation
// calculates derivative
matsu_state matsu_calc_deriv(double in, double t1recip, double t2recip,
double c, double b, double g,
matsu_state state)
{
double posX1 = POSPART(state.x1);
double posX2 = POSPART(state.x2);
double posIn = POSPART(in);
double negIn = NEGPART(in);
matsu_state deriv;
deriv.x1 = (c - state.x1 - (b*(state.v1)) - (g*posX2) - posIn) * t1recip;
deriv.x2 = (c - state.x2 - (b*(state.v2)) - (g*posX1) - negIn) * t1recip;
deriv.v1 = (posX1 - state.v1) * t2recip;
deriv.v2 = (posX2 - state.v2) * t2recip;
return deriv;
}
//helper function for moving the state forward by derivative*step
matsu_state matsu_eulerStep(matsu_state init, matsu_state deriv, int step)
{
matsu_state newVal;
newVal.x1 = init.x1 + (deriv.x1*step);
newVal.x2 = init.x2 + (deriv.x2*step);
newVal.v1 = init.v1 + (deriv.v1*step);
newVal.v2 = init.v2 + (deriv.v2*step);
return newVal;
}
// single RK4 step
void matsu_rkStep (double in, double t1recip, double t2recip,
double c, double b, double g,
matsu_state *state, int step){
matsu_state k1, k2, k3, k4;
k1=matsu_calc_deriv(in, t1recip, t2recip, c, b, g,
*state);
k2=matsu_calc_deriv(in, t1recip, t2recip, c, b, g,
matsu_eulerStep(*state, k1,step*0.5));
k3=matsu_calc_deriv(in, t1recip, t2recip, c, b, g,
matsu_eulerStep(*state, k2,step*0.5));
k4=matsu_calc_deriv(in, t1recip, t2recip, c, b, g,
matsu_eulerStep(*state, k3,step));
state->x1 = state->x1 + (k1.x1 + 2*(k2.x1+k3.x1) +k4.x1)*(1.0/6.0)*step;
state->x2 = state->x2 + (k1.x2 + 2*(k2.x2+k3.x2) +k4.x2)*(1.0/6.0)*step;
state->v1 = state->v1 + (k1.v1 + 2*(k2.v1+k3.v1) +k4.v1)*(1.0/6.0)*step;
state->v2 = state->v2 + (k1.v2 + 2*(k2.v2+k3.v2) +k4.v2)*(1.0/6.0)*step;
}
// wrapper to loop Rk4 step enough times to move forward by 1
double matsu_calc_nextVal_RK(double in, double t1, double t2,
double c, double b, double g,
matsu_state *state)
{
double step = 1.0 / STEP_DIVISION;
double t1recip = 1.0/t1;
double t2recip = 1.0/t2;
for(int i=0; i<STEP_DIVISION; i++){
matsu_rkStep(in, t1recip, t2recip, c, b, g, state, step);
}
return POSPART(state->x1) - POSPART(state->x2);
}
There will be no reason why you should use int as type of step in arguments of matsu_eulerStep and matsu_rkStep. You should use double just like other parameters.