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:
omega = 94 * pi / 30 ;
do = [1 4 7 11 18] ;
N = length(do) - 1 ;
tmax = 50 ;
tspan = [0 tmax];
x0 = [zeros(1,N) repmat(omega, 1, N)] ;
options = odeset('Mass', [eye(N) zeros(N) ; zeros(N) J]) ;
[t, xSol] = ode45(@(t, th) odefcn(t, th, thx, K, Tg, N) , tspan, x0, options) ;
subplot(2, 1, 1)
subplot(2, 1, 2)
function dxdt = odefcn(~, th, thx, K, Tg, N)
m = N + 1 ;
n = N * 2 ;
dxdt = zeros(n, 1) ;
T = diag(interp1(thx, Tg', wrapToPi(th(1:N)))) ;
dxdt(1:N) = th(m:n) ;
dxdt(m:n) = - K * th(1:N) + T;
I have tried different solvers and I am getting the same unexpected result.
Thanks for your help in advance,