ode45 is running an infinite loop

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 

Answer 1

[Star Strider]

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 

.

[Nicole]

Do you know what would go on the y axis of the first graph? Would it be what the x values represent?

[Nicole]

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!

[Walter Roberson]

The x values of the plots are time.

The y axis of the first plot is whatever x0(1) represents.

NrSp = size(x,2);

makes NrSp into the number of columns of the returned value x -- which will be the same thing as numel(x0), the number of states.

The for loop is then plotting one state at a time.

[Star Strider]

‘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.)

.