ode45 is running an infinite loop

4 views (last 30 days)
Nicole
Nicole on 23 Nov 2022
Commented: Star Strider on 29 Nov 2022
My ode45 keeps running an infinite loop. My code is shown below.
tspan = [0, 10]
x0 = [0;0;0;0;0]
T1 = 360000
V1 = 70
[t x] = ode45(@(t,x) odefun(t,x,T1,V1),tspan,x0)
plot(t,x)
function dxdt = odefun (t,x,T1,V1)
u = [T1;V1]
m = 18130.59;
k_t = 80363.83655;
c_aero = 4561.755;
c_t = 4561.755;
k_b = 62.86;
R = 1.41;%resistance
L = 0.000644;
I1 = 874897.2;
I2 = 0.000027525;
N = 149;%gear ratio
e= 2.1e11;
d_l = 0.7736;
d_h = 0.007;
l_l = 3.5;
l_h=0.3;
area_l = pi*(d_l/2)^2;
area_h = pi*(d_h/2)^2;
k_l = (e*area_l)/l_l;
k_h = (e*area_h)/l_h;
A = [0 1 0 0 0
(-N*k_h)/(I1) -c_aero/(I1) k_h/(I1) 0 0
0 0 0 1 0
(N*k_h)/(I2) 0 -k_h/I2 -c_t/I2 k_t/I2
0 0 0 -k_b/L -R/L]
B = [0 0
1/(I1) 0
0 0
0 0
0 -1/L]
C = eye(5)
D = [0]
dxdt = A*x+B*u
end

Answers (1)

Star Strider
Star Strider on 23 Nov 2022
The system is ‘stiff’ bercause of the large differences in the parameter values.
Use ode15s to get results —
tspan = [0, 10];
x0 = [0;0;0;0;0];
T1 = 360000;
V1 = 70;
[t x] = ode15s(@(t,x) odefun(t,x,T1,V1),tspan,x0)
t = 1000×1
1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
x = 1000×5
0 0 0 0 0 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0001 0.0000 0.0000 -0.0000 -0.0001 -0.0001 0.0000 0.0000 -0.0000 -0.0002 -0.0001 0.0000 0.0000 -0.0000 -0.0004 -0.0002 0.0000 0.0000 -0.0000 -0.0007 -0.0002 0.0000 0.0000 -0.0000 -0.0011 -0.0003
figure
plot(t,x)
grid
xlabel('Time')
NrSp = size(x,2);
figure
for k = 1:NrSp
subplot(NrSp,1,k)
plot(t, x(:,k))
grid
title(sprintf('x_%d',k))
end
xlabel('Time')
function dxdt = odefun (t,x,T1,V1)
u = [T1;V1];
m = 18130.59;
k_t = 80363.83655;
c_aero = 4561.755;
c_t = 4561.755;
k_b = 62.86;
R = 1.41;%resistance
L = 0.000644;
I1 = 874897.2;
I2 = 0.000027525;
N = 149;%gear ratio
e= 2.1e11;
d_l = 0.7736;
d_h = 0.007;
l_l = 3.5;
l_h=0.3;
area_l = pi*(d_l/2)^2;
area_h = pi*(d_h/2)^2;
k_l = (e*area_l)/l_l;
k_h = (e*area_h)/l_h;
A = [0 1 0 0 0
(-N*k_h)/(I1) -c_aero/(I1) k_h/(I1) 0 0
0 0 0 1 0
(N*k_h)/(I2) 0 -k_h/I2 -c_t/I2 k_t/I2
0 0 0 -k_b/L -R/L];
B = [0 0
1/(I1) 0
0 0
0 0
0 -1/L];
C = eye(5);
D = [0];
dxdt = A*x+B*u;
end
.
  4 Comments
Star Strider
Star Strider on 29 Nov 2022
Do you know what would go on the y axis of the first graph? Would it be what the x values represent?
I am not certain what the ‘x’ values represent, or what their units are, so I left the ylabel blank in the combined plot. If some of them asre positions, others velocities, and still others accelerations, it would be difficult to describe all of them in one axis label. (Generically, ‘Amplitude’ could be appropriate.) They could all be labelled specifically in the subplot array.
Would it be possible if you could explain the for loop and what it is doing and what NrSp is? It would be greatly appreciated!
Sure! The for loop creates an array of subplots, plotting each column of ‘x’ in a different subplot. (That is where labeling the y-axes would be appropriate, along with the units. They could also have individual titles and the subplot array could have one sgtitle.)
The ‘NrSp’ variable is the number of subplots, equal to the column size of ‘x’. (I chose subplot here to display the individual variables. Other options are stackedplot and tiledlayout. Use whatever best suits your needs. The syntax for them differs from the syntax for subplot, so my code would have to be changed slightly to use them.)
.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by