MATLAB: Using ODE45 to solve simultaneously independent ODE’s vs. calling the function twice – different results…

MATLABode45

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!

Best Answer

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.
Related Question