using ODE45 to solve simultaneously independent ODE's vs. calling the function twice - different results...

Hi,

i am trying to solve for the displacement of a spring & mass with rotating eccentric mass in X and Y directions.

If i solve it using the following code:

[t,y] = ode45(@odefun, tspan, ic);

    function dydt = odefun(t,y)
        %SDOF (Uncoupled) Forced (Unbalance), Undamped

        dydt = zeros(4,1);
        dydt(1) = y(2);
        dydt(2) = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);
        dydt(3) = y(4);
        dydt(4) = -(k/M)*y(3)+(m/M)*e*w^2*sin(w*t);

    end

I get one result (which seems incorrect for the second equation).

If i solve it using the following code (which effectively separates it into two independent ODEs:

[t,y] = ode45(@odefun, tspan, ic);

    function dydt = odefun(t,y)
        %SDOF (Uncoupled) Forced (Unbalance), Undamped

        dydt = zeros(2,1);
        dydt(1) = y(2);
        dydt(2) = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);

    end

[t2,y2] = ode45(@odefun2, tspan, ic2);

    function dydt2 = odefun2(t,y2)

        dydt2 = zeros(2,1);
        dydt2(1) = y2(2);
        dydt2(2) = -(k/M)*y2(1)+(m/M)*e*w^2*sin(w*t);

    end

I get a different result, which seems more in line with what i expect.

Why would i get different results for these two methods?

See attached images..

Any help would be greatly appreciated!

[James Tursa]

What initial conditions are you using? When I run the code with the ic from your other post and M=1 I get plots that are very nearly right on top of each other (with the exception that the output sample times are different).

[Tyler]

k = 104000 M = 200 m = 1/16 w = 1 e = 1

InitialX,InitialVx,InitialY,InitialVy = 0

Answer 1

[Jan]

This is a problem of the accuracy only. Decrease the tolerance to get equivalent trajectories for the simultaneous and separate evaluation:

function YourFcn

for method = 1:2
  if method == 1
    opt = odeset('RelTol', 1e-3, 'AbsTol', 1e-6);
  else
    opt = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
  end

    % Solve simultaneously:
    tspan = [0, 10];
    ic    = [0,0,0,0];
    [t,y] = ode45(@odefun, tspan, ic, opt);

    figure;
    for k = 1:4
      subplot(2,4,k);
      plot(t, y(:,k));
    end

    % Solve separately:
    ic    = [0,0];
    [t,y] = ode45(@odefun1, tspan, ic, opt);
    subplot(2,4,5);
    plot(t, y(:,1));
    subplot(2,4,6);
    plot(t, y(:,2));

    [t,y] = ode45(@odefun2, tspan, ic, opt);
    subplot(2,4,7);
    plot(t, y(:,1));
    subplot(2,4,8);
    plot(t, y(:,2));
  end
  end

function dydt = odefun(t,y)
k = 104000;
M = 200;
m = 1/16;
w = 1;
e = 1;
dydt = [y(2); ...
        -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t); ...
        y(4); ...
        -(k/M)*y(3)+(m/M)*e*w^2*sin(w*t)];
end

function dydt = odefun1(t,y)
k = 104000;
M = 200;
m = 1/16;
w = 1;
e = 1;
dydt = [y(2); ...
        -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t)];
end

function dydt = odefun2(t,y)
k = 104000;
M = 200;
m = 1/16;
w = 1;
e = 1;
dydt = [y(2); ...
        -(k/M)*y(1)+(m/M)*e*w^2*sin(w*t)];
end

High tolerances - different results: First row is the simultaneous solution, second row the separated:

Lower tolerances, which means a more accurate integration - equivalent results:

With the lower tolerances the trajectory shows the expected periodic behavior also, while for the less accurate integration the local accumulated discreatization error dominates the solutions. It is simply too coarse.

Note that the default of the absolute tolerance of 1e-6 is in the magnitude of your values. Then this is not a meaningful limit to determine a proper stepsize.

[Tyler]

Excellent! Thank you so much for your good work!

[Jan]

You are welcome.

[Chenko DUSHIMIMANA]

Hello, Thanks the detail explanation about ode45. What if the differential equations to be solved by ode45 are interconnected(i.e. they depend on each other, that is, to find the solution of one system of equations(say, system 1) requires the knowledge of the solution from the other system of equations(say, system 2)). Where, the "tspan" for both system may differ or be similar. How Can we solve these two systems? simultaneously or separately?

[Jan]

For interdependencies you have to solve the equations together. Then there is the special case, that the two systems have a very different frequency behaviour, e.g. two coupled pendulums and the eigenfrequency of one is 1e8 times higher than the other. In this case, the step size control of ODE45 considers the high frequency part of the system only and the accumulated rounding errors can dominate the solution for the low frequency part. There is no smart solution to handle this problem, because the Runge-Kutta-one step solvers are not designed for this case. You need a stiff solver for a reliable solution.

This is a problem for the simulation of a car: If you simulate the body, the wheels and consider the engine as a black box producing force, you will not catch the high-frequency vibrations which cause discomfort and noise. A complete model including all small parts is very expensive and you have to choose the integrator carefully. Even then the parameters of the integrator can influence the results massively.

I've seen (too) many scientific calculations, in which a complex system of differential equation was "solved" by a standard integrator like ODE45. Without an analysis of the stiffness (e.g. caused be extremely different eigenfrequencies) and the sensitivity (small variation of initial values and parameters), the result is not reliable or valid. Numerical maths is a profession and it takes some years to learn how to apply and control the different methods. So if you have to integrate a conmplicated system, ask a scientist working in the field of numerical maths.

Btw. it was a development team of VW, who ran a simulation of a detailed car model using an ODE45 integrator written by their own. They failed tremendously due to bugs in the stepsize control and the stiffness of the system. Fortunately this was in the 1990'th and today experienced mathematicians are members of the team.

Answer 2

[Abraham Boayue]

Answers have been inserted into your codes (I just modified them a bit).

function Plot_odes 
T = 10;
N = 1000;
t = 0 :T/(N-1):T;
ic = [0 0]; 
[t,x]= ode45(@odefun, t,ic );
[t,y]= ode45(@odefun2, t,ic );

figure
plot(t ,x(:,1),'linewidth',2,'color','r')
hold on
plot(t ,y(:,1),'linewidth',2,'color','b')
legend('x','y');
a = title('x(t) y(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 T]');
set(a,'Fontsize',14);
xlim([0 T])
grid

% Your two ODEs are completely independent of each other and
% therefore, the solution of one does not interfere with the other. In other
% words, they are not simultineous equations; you actually tricked ode45
% into believing that it was solving simultaneous equations in the first
% case. However, the two functions (odefun and odefun2) below are valid solutions of
% the odes because the vectors dxdt = [x' x''], dydt = [y' y''] and their
% solution vectors, X = [ x1 x2] and Y = [y1 y2] (where x1 = x(t) as the
% required solution for the displacement in the x direction and x2 its
% derivative), are uniform vectors with each element dependent on the other. On the otherhand, in your first method,
% you had dydt = [x' x'' y' y''] and Y = [x1 x2 y1 y2] which can not be right 
% because x and y are independent variables.

 function dxdt = odefun(t,y)
          %SDOF (Uncoupled) Forced (Unbalance), Undamped
          k = 104000;
          M = 200; 
          m = 1/16; 
          w = 1; 
          e = 1;
          dxdt_1 = y(2);
          dxdt_2 = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);
          dxdt =[dxdt_1; dxdt_2]; 

          function dydt = odefun2(t,y2)
          k = 104000;
          M = 200; 
          m = 1/16; 
          w = 1; 
          e = 1;
          dydt_1 = y2(2);
          dydt_2 = -(k/M)*y2(1)+(m/M)*e*w^2*sin(w*t);
          dydt =[dydt_1; dydt_2];

Here is your first method.

function Plot_odes 
T = 10;
N = 1000;
t = 0 :T/(N-1):T;
ic = [0 0 0 0]; % Initial conditions 
[t,z]= ode45(@odefun, t,ic );

figure
plot(t ,z(:,1),'linewidth',2,'color','r')
hold on
plot(t ,z(:,3),'linewidth',2,'color','b')
legend('x','y');
a = title('x(t) y(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 7]');
set(a,'Fontsize',14);
xlim([0 T])
grid

% Even though you were able to get some results, however, 
% they are misleading because ode45 sees
% dydt = [y' y'' y''' y''''] and gives Y = [y1 y2 y3 y4].
% See above for further explanations
      function dydt = odefun(t,y)
          %SDOF (Uncoupled) Forced (Unbalance), Undamped
          k = 104000;
          M = 200; 
          m = 1/16; 
          w = 1; 
          e = 1;
          dydt_1 = y(2);
          dydt_2 = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);
          dydt_3 = y(4);
          dydt_4 = -(k/M)*y(3)+(m/M)*e*w^2*sin(w*t);
          dydt =[dydt_1; dydt_2;dydt_3; dydt_4];

[Jan]

A Runge-Kutta integrator like ODE45 does not draw the conclusion that the 3rd component is the derivative of the 2nd one. The relation between the components is defined in the equation to be integrated only. It is valid to integrate independent components, as long as the stepsize control gets sufficient parameters for the tolerances. The only drawback is that the stepsize controller cannot run the integration with the optimal step width for the simultaneous solution in opposite to the separated functions.

[Abraham Boayue]

Great job Jan! You have just thought me a valuable lesson.

[Chenko DUSHIMIMANA]

Hello, Thanks the detail explanation about ode45. What if the differential equations to be solved by ode45 are interconnected(i.e. they depend on each other, that is, to find the solution of one system of equations(say, system 1) requires the knowledge of the solution from the other system of equations(say, system 2)). Where, the "tspan" for both system may differ or be similar. How Can we solve these two systems? simultaneously or separately?