MATLAB: Ode45 unexpected behaviour for initial conditions = 0

differential equationsinitial conditionsode23ode23tode45

Hi all,
I have a code where I am using it to solve the following differential equation:
In specific, the above is a 4×4 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.

Best Answer

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