ode45 unexpected behaviour for initial conditions = 0

Hi all,

I have a code where I am using it to solve the following differential equation:

In specific, the above is a 4x4 system of equations where J is the polar moment of inertia (a diagonal matrix), K is a stiffness matrix, and T_g(th) is the excitation torque on the system.

What I want to do is simply solve the system for some initial conditions that I set myself. However, just as a sanity check, I solve the equation with all initial conditions set to zero. Naturally, this should give a zero response for the equation's solution, however it doesn't which I find quite odd. Hence with that result, I cannot trust that I will have a correct answer with any initial condition.

So any ideas of what might be going wrong?

Attached you can find the .mat file that I draw all my inputs from and my ode45 code is below:

clear
clc

% Input data

    % Input Engine Data
    omega = 94 * pi / 30 ;                  % Rotational Speed as initial condition

    % Degrees of Freedom
    do = [1 4 7 11 18] ;                    % DOF Split vector
    N = length(do) - 1 ;                    % Number of DOFs

    % ODE Solution
    tmax = 50 ;                             % Solution time (s)
    
    % Load ODE Matrices
    load('Struct.mat')
          
%% ODE Solution

% ODE Initial Conditions
    % Time
    tspan = [0 tmax];
    
    % Initial conditions
    x0 = [zeros(1,N) repmat(omega, 1, N)] ; 
    
    % ODE options
    options = odeset('Mass', [eye(N) zeros(N) ; zeros(N) J]) ;

    % Solution
    [t, xSol] = ode45(@(t, th) odefcn(t, th, thx, K, Tg, N) , tspan, x0, options) ;
 %%
    figure
    subplot(2, 1, 1)
    plot(t, xSol(:,1:N))
    title('Displacement')
    subplot(2, 1, 2) 
    plot(t, xSol(:,N+1:end))
    title('Velocity')
   
%% ODE
    
function dxdt = odefcn(~, th, thx, K, Tg, N) 

% Indices
m = N + 1 ;
n = N * 2 ;

% Preallocate
dxdt = zeros(n, 1) ;

% Torques
T = diag(interp1(thx, Tg', wrapToPi(th(1:N)))) ;    % Gas Torque

% ODE Equation
dxdt(1:N) = th(m:n) ;
dxdt(m:n) = - K * th(1:N) + T; 

end
        
    
    

P.S.

I have tried different solvers and I am getting the same unexpected result.

Thanks for your help in advance,

KMT.

Answer 1

[Jon]

I think you have a basic conceptual error in your ode definition.

You should have an overall state vector x = [x1;x2] whre x1 = theta, x2 = d(theta)/dt so your ode's look like

dx1/dt = x2

dx2/dt = -K*theta + T(x)

you have instead dx1/dt = theta, which is not correct

[KostasK]

Hi Jon,

Thanks for your response. If I understand correctly what you mean, I should have done this part correctly. The reason being the variable that appears in the code as th(m:n) is d(theta)/dt, and the variable that appears in the code as th(1:N) is theta.

The reason that I am saying this is because I have tested it on a simple damped mass-spring system and it behaves correctly (attached you can find the code for that if you are interested).

I did spend some time playing with this variable in the code, so in case th(m:n) shoudln't be there cause its wrong, I am not aware with that else to substitute it with in the code.

[Jon]

Yes, I see that it looks like your equations are equivalent to mine. I didn't look carefully enough at your indices. For the future, the code might be a little more obvious if you did something like

function dxdt = odefcn(~, th, thx, K, Tg, N) 


% Torques
T = diag(interp1(thx, Tg', wrapToPi(th(1:N)))) ;    % Gas Torque
% ODE Equation

% assign local variables for readability
theta = th(1:N);
omega = th(N+1:end);

% compute derivatives
dthetadt = omega;
domegadt = -K*theta + T;

% make overall vector of derivatives
dxdt = [dthetadt;domegadt];

When you do your test setting all of the initial conditions to zero, do you somehow also force the torque to zero. Otherwise, if I understand your system correctly, wouldn't a non-zero torque cause the system to accelerate even if the initial conditions were zero?

[KostasK]

''...I understand your system correctly, wouldn't a non-zero torque cause the system to accelerate even if the initial conditions were zero?''

This is a valid point, and I have been thinking of that as well, so here is my logic: The torque is strictly a function of theta and not a function of time (which is the independent variable). In addition, the torque function is such that T(0)=0. In which case I figure that if theta begins as zero, there will be no reason for the excitation torque to become anything different than zero as well. Hence the response should be zero as well. By the same token, if the excitation torque was time dependent, then yes as you said the response would be non-zero as well.

[KostasK]

Indeed you are correct however! I just had a look and as the function is approximated near zero, it is not exactly equal to zero when it should. Given how small the aplitudes of the oscillation are I would say that this is what could be causing this