Why does my ode45 function not run?
1 view (last 30 days)
Show older comments
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 Comments
Answers (1)
Sam Chak
on 7 Apr 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 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!