Why does my ode45 function not run?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (1)
Sam Chak
le 7 Avr 2023
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
0 commentaires
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!