Using ODE1 in Simulink yields different results for an equivalent model run in a for loop

I've run into an interesting situation with Simulink. I am using Simulink to evaluate alternative numerical integration routines for some time-varying ordinary differential equations. The Simulink model that I have built is provided in the attached screen shot. The block, quaterionKinematics, encodes the continuous-time, time-varying differential equation that governs the evolution of quaternion kinematics. As you can see, the output of this block is the time derivative of the quaternions, we integrate to get q(t), and then perform a normalization step.

function q_dot = quaternionKinematics(wx, wy, wz, q)
 
q_dot = 0.5*[0, -wx, -wy, -wz;...
            wx, 0, wz, -wy;...
            wy, -wz, 0, wx;...
            wz, wy, -wx, 0]*q;
 
end

The gyroscope data, given by the wx, wy, wz is loaded from my workspace. In evaluating these results against alternative integraiton schemes that I will have to implement in another programming language, I have built the equivalent for loops that perform this same integration. Equivalent for loop:

%% Euler Integration

h = Ts; %sample period
Omega(:,:,1) = 0.5*[0 -wx(1) -wy(1) -wz(1);...
        wx(1) 0 wz(1) -wy(1);...
        wy(1) -wz(1) 0 wx(1);...
        wz(1) wy(1) -wx(1) 0];
 
q_last = initQuat;    
Omega_last = Omega(:,:,1);
for i = 2:length(time_i)
    
    Omega(:,:,i) = 0.5*[0 -wx(i) -wy(i) -wz(i);...
                        wx(i) 0 wz(i) -wy(i);...
                        wy(i) -wz(i) 0 wx(i);...
                        wz(i) wy(i) -wx(i) 0];
    
    q_Euler(:,i) = q_last + h*Omega_last*q_last;
    
    q_Euler(:,i) = q_Euler(:,i)/norm(q_Euler(:,i),2);
    
    q_last = q_Euler(:,i);
    Omega_last = Omega(:,:,i);
    
    
end

When using ODE1, which is the forward Euler method that I have depicted above, my results disagree around certain points in my data...and by quite a bit. I have included a figure that shows the difference between one of the propagated quaternion elements.

When I move toward higher-order integration schemes, such as ODE2 (Heun), the output of Simulink tends to track the output of what should be an equivalent for loop a bit better, however, there is still inconsitency between the two outputs, which I would have thought should be equivalent. The most significant disagreement seems to be specific to the results that I am getting with ODE1.

Can you provide some guidance on why this may be? I am hesitant to use the results that I am obtaining from Simulink in making design decisions, now. Could this be related to how I am encoding my differential equations, where I have used an embedded Matlab function? This approach has not cuased problems for me in the past, and again, gross disagreement in results really only seem to be specific to use of ODE1 as my solver.

[James Tursa]

Is the Simulink version internally using a variable step size?

[Chris D'Angelo]

Since ODE1 is a fixed step size solver, I would hope not. I could see issues arising if - internally - Simulink did not properly handle the sample period that is provided.

[Chris D'Angelo]

Something that I have not tried...since one of my primary interests here was in evaluating alternative integration schemes, was the manual codification of the discrete-time ODE: that would defeat the purpose, however, of relying on Simulink to transform my system equations into these alternative, higher-order forms. Hence my decision to encode the continuous-time version of my differential equations!

Answer 1

[Paul]

I think some of the equations are out of order. I got a match (obviously using my inputs) with this:

h = Ts; %sample period
q_Euler = zeros(4,numel(out.tout));
q_Eulerdot = zeros(4,numel(out.tout));
q_Euler(:,1) = out.qsim(1,:).';  % out.qsim is the normalized quaternion, output from simulink model
q_int = out.qsim1(1,:).';  % out.qsim1 is the output of the integrator (non-normalized quaternion) output from simulink model

for i = 1:length(out.tout)
    
    q_Euler(:,i) = q_int/norm(q_int);  % normalize the integrator output
    
    % compute Omega with wx, wy wz outputs from Simulink model
    Omega = 0.5*[0 -out.wx(i) -out.wy(i) -out.wz(i);...
        out.wx(i) 0 out.wz(i) -out.wy(i);...
        out.wy(i) -out.wz(i) 0 out.wx(i);...
        out.wz(i) out.wy(i) -out.wx(i) 0];
    
    q_Eulerdot(:,i) = Omega*q_Euler(:,i);

    q_int = q_int + h*q_Eulerdot(:,i);  % propagate integrator to next time step
    
end

It will be important that the quat_normalization function implements the /norm() as in the code above.

Historical note: a long time ago I ran into a similar problem where Simulink didn't exactly match the Matlab implementation, and was told by Tech Support that Simulink and Matlab could be linked with different libraries, so innoccuous operations like matrix multiply might not even give the exact same result. Don't know if that's still true.