Accuracy problems

[mohit-s]

3 sec simulation of RLC with odeint

3 sec simulation of RLC with matlab

0.1 sec simulation of RL with odeint

0.1 sec simulation of RL with matlab

My original model is more complex, but to debug the problem, I am using this simplest system: const double M = 9.4e-3; const double cm = 0.12; const double f = 24.7;

typedef boost::array<double, 2> state_type;

void rlc (const state_type &x, state_type &dxdt, double t) { dxdt[0] = x[1]; dxdt[1] = -(cm/M)*x[1] + 1; }

void write_rlc (const state_type &x, const double t) { cout << t << '\t' << x[0] << '\t' << endl; }

int main(int argc, char **argv) { state_type xnow = {0.0, 0.0}; double tcur = 0.0; double tstep = 0.001/f; double tfinal = 0.1;

runge_kutta4 stepperRK4; integrate_const(stepperRK4, rlc, xnow, tcur, tfinal, 0.1*tstep, write_rlc); }

The result does not match with MATLAB solution: function ydot = piezo_rlc(t,y)

kp = 226; cm = 1.7*0.12; TH = 0.001914; g = 0.98; M = 9.4e-3; C = 39.2e-9; C_st = 1e-6; %f=29.2; f=24.7;

%ydot = [y(2); -(cm/M)_y(2) - (kp/M)_y(1) + g_sin(2_pi_f_t)]; yo(1) = y(2); yo(2) = -(cm/M)_y(2) + 1; %yo(2) = -(cm/M)_y(2) - (kp/M)_y(1) + g_sin(2_pi_f*t); ydot = [yo(1); yo(2)];

end

[to,yo] = ode45(@piezo_rlc, [0, 0.1], [0 0]);

Note that when I use the full rlc electrical system, the frequency is also different from MATLAB solution. I have also used integrate_adaptive, which gives the same result as integrate_const.

[mariomulansky]

the matlab code uses a runge kutta method with error control (ode45), while in C++ you use the basic rk4 which doesnt have error control. if you switch to something else like dopri5 in the odeint example the results might get more similar.

[mohit-s]

Hi, thanks for the quick reply. I modified the code to try both error control methods. The results are as before. I have also used controlled_stepper. Also ode23 in matlab gives correct result.

typedef runge_kutta_cash_karp54< state_type > error_stepper_type_CK54;
error_stepper_type_CK54 stepperCK54;
integrate_const( stepperCK54, rlc, xnow, tcur, tfinal, 0.1*tstep, write_rlc);

typedef runge_kutta_dopri5< state_type > error_stepper_type_D5;
error_stepper_type_D5 stepperD5;
integrate_const( stepperD5, rlc, xnow, tcur, tfinal, 0.1*tstep, write_rlc);

[Note the graph is still the same as original post. matlab graph grows slower and has more bend.]

[headmyshoulder]

The constants for cm are different for matlab and c++. Could this explain the difference?

On 29.06.2014 21:56, mohit-s wrote:

3 sec simulation of RLC with odeint odeint_rlc_x1 https://cloud.githubusercontent.com/assets/5552614/3424261/fdb25d8c-ffc6-11e3-9316-e296575a7120.png

3 sec simulation of RLC with matlab matlab_rlc_x1 https://cloud.githubusercontent.com/assets/5552614/3424262/00d386c6-ffc7-11e3-93b1-421ccf6e5938.png

0.1 sec simulation of RL with odeint odeint_x1 https://cloud.githubusercontent.com/assets/5552614/3424265/126e8674-ffc7-11e3-97c4-18cab337f3e6.png

0.1 sec simulation of RL with matlab matlab_x1 https://cloud.githubusercontent.com/assets/5552614/3424266/15a4b250-ffc7-11e3-9e31-2d86cd3f31bc.png

My original model is more complex, but to debug the problem, I am using this simplest system: const double M = 9.4e-3; const double cm = 0.12; const double f = 24.7;

typedef boost::array state_type;

void rlc (const state_type &x, state_type &dxdt, double t) { dxdt[0] = x[1]; dxdt[1] = -(cm/M)*x[1] + 1; }

void write_rlc (const state_type &x, const double t) { cout << t << '\t' << x[0] << '\t' << endl; }

int main(int argc, char **argv) { state_type xnow = {0.0, 0.0}; double tcur = 0.0; double tstep = 0.001/f; double tfinal = 0.1;

runge_kutta4 stepperRK4; integrate_const(stepperRK4, rlc, xnow, tcur, tfinal, 0.1*tstep, write_rlc); }

The result does not match with MATLAB solution: function ydot = piezo_rlc(t,y)

kp = 226; cm = 1.7*0.12; TH = 0.001914; g = 0.98; M = 9.4e-3; C = 39.2e-9; C_st = 1e-6; %f=29.2; f=24.7;

%ydot = [y(2); -(cm/M)_y(2) - (kp/M)_y(1) + g_sin(2_pi_f_t)]; yo(1) = y(2); yo(2) = -(cm/M)_y(2) + 1; %yo(2) = -(cm/M)_y(2) - (kp/M)_y(1) + g_sin(2_pi_f*t); ydot = [yo(1); yo(2)];

end

    [to,yo] = ode45(@piezo_rlc, [0, 0.1], [0 0]);

Note that when I use the full rlc electrical system, the frequency is also different from MATLAB solution. I have also used integrate_adaptive, which gives the same result as integrate_const.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/132.

[mohit-s]

You are right. I got the constants mixed up. The more complex model also matches. Thanks for your time, and sorry for the distraction.