Why does my ode45 function not run?

This code just keeps loading. It does not run. If I change my initial values to [0; 0; 0; 1;0;-1;0] it runs fine. But with any other intial values it just keeps loading. How do I fix this?

t = 0:0.01:(2*pi*4);

options=odeset('RelTol',1e-12,'AbsTol',1e-12);

[t,q] = ode45(@euler_eqns, t, [0.5; 0; 0; 0.866;0;-1;0], options);

K = q(:,1).^2 + q(:,2).^2 + q(:,3).^2 + q(:,4).^2;

K0 = 1;

c = K - K0;

plot(t,K);

xlabel('Time (s)');

ylabel('K');

plot(t,c);

xlabel('Time (s)');

ylabel('K-K0');

function dqwdt = euler_eqns(t,qwd, q)

q = qwd(1:4);

w = qwd(5:7);

I = 500;

J = 125;

K = (I-J)/I;

T = 35;

ohm = 1;

qstar = 1/ohm;

dqdt = zeros(4,1);

dqdt(1) = pi*(q(4)*w(1) + q(3)*(qstar - w(2) - ohm) + q(2)*w(3));

dqdt(2) = pi*(q(3)*w(1) + q(4)*(w(2) - ohm - qstar) - q(1)*w(3));

dqdt(3) = pi*(-q(2)*w(1) - q(1)*(qstar - w(2) - ohm) + q(4)*w(3));

dqdt(4) = pi*(-q(1)*w(1) - q(2)*(w(2) - ohm - qstar) - q(3)*w(3));

dwdt = zeros(3,1);

dwdt(1) = 12*pi*K*(q(2)*q(3) + q(1)*q(4))*(1 - 2*q(1)^2 - 2*q(2)^2) ...

- 2*pi*K*w(2)*w(3) + 2*pi*qstar*w(3);

dwdt(2) = 0;

dwdt(3) = -24*pi*K*(q(3)*q(1) - q(2)*q(4))*(q(2)*q(3) + q(1)*q(4)) ...

+ 2*pi*K*w(1)*w(2) + 2*pi*qstar*w(1);

% Combine the time derivatives into a single vector

dqwdt = [dqdt; dwdt];

end

Answer 1

[Sam Chak]

Hi @Sneh

That's because the response of one of the states, exploded. This it cannot run to the end sec.

% t = 0:0.01:(2*pi*4);
t = 0:0.1:3.9;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[t, q]  = ode45(@euler_eqns, t, [0.5; 0; 0; 0.866; 0; -1; 0]);
% [t, q]  = ode45(@euler_eqns, t, [0.5; 0; 0; 0.866; 0; -1; 0], options);
plot(t, q(:,5))

% K  = q(:,1).^2 + q(:,2).^2 + q(:,3).^2 + q(:,4).^2;
% K0 = 1;
% c  = K - K0;
% plot(t, K);
% xlabel('Time (s)');
% ylabel('K');
% plot(t, c);
% xlabel('Time (s)');
% ylabel('K-K0');

function dqwdt = euler_eqns(t,qwd, q)
    q       = qwd(1:4);
    w       = qwd(5:7);
    
    I       = 500;
    J       = 125;
    K       = (I-J)/I;
    T       = 35;
    ohm     = 1;
    qstar   = 1/ohm;
    
    dqdt    = zeros(4,1);
    dqdt(1) = pi*( q(4)*w(1) + q(3)*(qstar - w(2) - ohm)   + q(2)*w(3));
    dqdt(2) = pi*( q(3)*w(1) + q(4)*(w(2)  - ohm  - qstar) - q(1)*w(3));
    dqdt(3) = pi*(-q(2)*w(1) - q(1)*(qstar - w(2) - ohm)   + q(4)*w(3));
    dqdt(4) = pi*(-q(1)*w(1) - q(2)*(w(2)  - ohm  - qstar) - q(3)*w(3));
    
    dwdt    = zeros(3,1);
    dwdt(1) =  12*pi*K*(q(2)*q(3) + q(1)*q(4))*(1 - 2*q(1)^2 - 2*q(2)^2) - 2*pi*K*w(2)*w(3) + 2*pi*qstar*w(3);
    dwdt(2) = 0;
    dwdt(3) = -24*pi*K*(q(3)*q(1) - q(2)*q(4))*(q(2)*q(3) + q(1)*q(4)) + 2*pi*K*w(1)*w(2) + 2*pi*qstar*w(1);
    
    % Combine the time derivatives into a single vector
    dqwdt = [dqdt; 
             dwdt];
    
end