Azzera filtri
Azzera filtri

Why does my ode45 function not run?

1 visualizzazione (ultimi 30 giorni)
Sneh
Sneh il 7 Apr 2023
Risposto: Sam Chak il 7 Apr 2023
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

Risposte (1)

Sam Chak
Sam Chak il 7 Apr 2023
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

Categorie

Scopri di più su Language Fundamentals in Help Center e File Exchange

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by