How many iterations does a ODE solver like ode45 do? What do the dimension of outputs mean?

Say I have a state space function

function  dx = sys(t, x)
   A = [0 1; -2 3]; B = [0;1]; K = [-1 -1]; 
   u = K*x;
   dx = A*x + B*u;
end

Which I will call with the following script:

tspan = [0 10];
iniCon = [1;-1]
[t, x] = ode45(@sys, tspan, iniCon);
plot(t,x)

This way, I have a two-variable function of x with an initial condition of [1;-1] which I want plotted over 0 to 10 time units.

When running this program, I get t and x with dimension 85x1 and 85x2, respectively.

My questions which are made explicit in its subquestions related to this example are:

1) Why do I get 85 rows in x?

1a)- Where is the number 85 coming from?

1b)- Why are rows and columns reverted? x started with 2 columns, now it is 2 rows.

2) How many iterations does ode45 do and why this amount?

2a) Similarly to (1a): From a timespan of [1 10] I get 85 timestamps, what is the relation and why?

Thanks in advance. See a plot below of plot(t,x).

Answer 1

[Walter Roberson]

When you pass in tspan of length 2, ode45() decides for itself when to output.

ode45() takes variable step lengths. At each step it evaluates the function at carefully chosen locations and use the results to make a prediction about where the function is going. It then takes the prediction and compares to the actual results. If the two are close enough then matlab "accepts" the point and records it internally, and increases the step size and goes back to make another prediction. If the two were not close enough then ode45 rejects the step and reduces the step size and tries again. It might fail several times in a row in places that the function is changing rapidly.

When the tspan is a vector of length 2 then each time it accepts a point, then by default it uses the predictions to generate 3 intermediate points between the previous and new points, and predicts the values there and emits those predictions and the accepted point.

So in your case about 28-ish points were accepted, 3 additional predictions for each one, plus the start point, gives 85 outputs.

The reason that the rows become columns on output is that plot() plots each column as a separate line, so it makes more sense to have each the different variables be different columns to make plot easier.

ode45 does not do iterations as such. It just keeps evolving the function and checking if the prediction is reasonable. It rejects if the function changes more rapidly than the prediction function.

One way of thinking of it is that it is making a series of quintic (order 5) polynomial approximations of the function. It uses as many tries as is needed to do that.

[DB]

Thank you, this makes it more clear. Why I'm asking this question is because my own script is also a state space system with x in 6x1, A 6x6, B 6x3, with u = k where k is a 3x1 vector calculated with some programming where each row consists of a very long polynomial. When I ran that script it took endless amounts of time to run to the point I stopped it. Would you have any tips or trick in running such a expansive system? I think I can already take that setting tspan to a length of 1 greatly reduces running time. I solely need a plot to see whether u=k asymptotically converges the system or not, so a few iterations should be sufficient.

[DB]

I see now that tspan must be length of 2, as length of 1 gives an error with the function handle, and tspan = [0 10 20] gives strange results.

[Bruno Luong]

The number of steps mainly depends on how stiff your system is

You can adjust the timestep rejection sensitiviyu by providing RealTol AbsTol, see https://uk.mathworks.com/help/matlab/math/summary-of-ode-options.html

You migh also change to stiff solver such as ode15s instead of ode45.

If you want to learn more about the rejection, read the paper in the bottom of ode45 doc page. You can edit ode45.mand see how the error is estimated:

fE = f1*e1 + f3*e3 + f4*e4 + f5*e5 + f6*e6 + f7*e7;

where f1, ..., f7 are various ode intermediate evaluations, e1, ..., e7 are some constant to estimate the error (sum to 0). If the error is larger than some quantity depending on RealTol AbsTol, the timestep is rejected and divided by half, ... until success.

[Walter Roberson]

tspan must have at least two elements. If it has 3 or more then each time it accepts a point at or after one of the times in the vector it makes a prediction about what the value would be at that time and emits it.

[Walter Roberson]

I suggest you switch to ode23s instead of ode45. You probably have a stiff system.

[Sam Chak]

Hi @DB

If you choose a different set of values for K, then the system won't blow up. For example, .

tspan  = [0 10];
iniCon = [1;-1];
[t, x] = ode45(@sys, tspan, iniCon);
plot(t, x)

A  = [0 1; -2 3]; 
B  = [0; 1]; 
K  = [1 -5]; 
AA = A + B*K    % A*x + B*(K*x) = (A + B*K)*x

AA = 2×2

0 1 -1 -2

eig(AA)

ans = 2×1

-1 -1

function  dx = sys(t, x)
   A  = [0 1; -2 3]; 
   B  = [0; 1]; 
   K  = [1 -5]; 
   u  = K*x;
   dx = A*x + B*u;
end