Is it possible to extract the time and two output variables from ODE45?

Hello,

I have written a function which holds state space equations for ODE45 to solve, being in the nature of:

function zdot = odess(t,z)

In my mainscript which calls upon this function and has the initial conditions for the ode, I am recalling the outputs by:

[t,z] = ode45(@odess, t_step, z0, opts);

However, is it possible to get a second variable from the 2nd code listed above? For example:

[t,z,new_output] = ode45(@odess, t_step, z0, opts);

Thanks for your time.

[Jan]

Where does the "new_output" come from? What do you expect in this variable?

[Tyler Bikaun]

Sorry Jan, the syntax in the 3rd code above is just off of my head and not what I actually expect to be able to do. I want to be able to extract z(1) from my odescript (see comments below), I thought that that may be a way to achieve it.

Answer 1

[Torsten]

The easiest way is to reculculate new_output from t and z after the call

[t,z] = ode45(@odess, t_step, z0, opts);

Best wishes

Torsten.

[Tyler Bikaun]

I'm not quite sure how I'd be able to recalculate the output I require from the outputs given by ODE45.

My code within 'odess' is:

function zdot = odess(t,z)
% defining global variables
global k_mb r_p r_g I_p I_g m_p m_g V

%%state space equations are defined below:
zdot(2) = ((r_p*k_mb)/I_p)*(-r_p*z(1)+r_g*z(3)+z(5)-z(7));  % Eq 1ss
zdot(6) = (k_mb/m_p)*(r_p*z(1)-r_g*z(3)-z(5)+z(7));         % Eq 2ss
zdot(4) = ((r_g*k_mb)/I_g)*(-r_g*z(3)+r_p*z(1)+z(7)-z(5));  % Eq 3ss
zdot(8) = (k_mb/m_g)*(-r_p*z(1)+r_g*z(3)-z(7)+z(5));        % Eq 4ss

zdot(1) = z(2); % Eq 5ss
zdot(3) = z(4); % Eq 6ss
zdot(5) = z(6); % Eq 7ss
zdot(7) = z(8); % Eq 8ss
% Outputs
zdot = zdot';

I am getting my time step and each of the zdot's outputted but I also want the state variable z(1) to be an output as well.

I thought of doing:

zdot(9) = z(1);

But I think it doesn't work that easily..

[Torsten]

But z(1) is automatically included in the matrix z you get after the call to ode45.

Try

[t,z] = ode45(@odess, t_step, z0, opts);
plot (t,z(:,1));

to get a plot of z(1) over time.

Best wishes

Torsten.

[Walter Roberson]

No, z(1) input to the function is not reliably included in the z output . The ode outputs are not necessarily for the same time points that the inputs were evaluated at: the ode routines make projections of values at time points. Also, if you supplied a list of times to evaluate at (rather than a vector with start and end time) then the outputs are only at the times you indicated but the times evaluated at are whatever is necessary to achieve the precision specified.

[Torsten]

z is the solution variable vector - and each solution variable (in this case z(1)) is included in the output at the times t if you call ode45 as

[t,z] = ode45(@odess, t_step, z0, opts);

Best wishes

Torsten.

[Jan]

@Torsten: I agree. The zdot are not written to the output, because they are used internally only for the integration. Then e.g. "zdot(1) = z(2)" does not matter here.

I still do not understand, what "new_output" should be. Therefore I cannot post an own answer and even cannot reconsider, what the otehr answers are talking of.

[Tyler Bikaun]

z(1) isn't included in the output matrix "z" from [t,z] = ode45(@odess,t_step,z0, opts), i think what is included are only zdots.

So, I've added zdot(9) = z(1) as a new state space equation within the odess function and updated z0 for an additional initial condition.

As z(1) is supposed to be angular displacement of a gear, the output seems to be reasonable as it is not oscillating in value like the others outputs (but I cannot be 100% sure of this, could be a coincidence). However, I'm not sure if the inclusion of that line of code actually is integrated by ode45 when its ran, or just outputs the value of z(1), if its the former, then its not outputting what is required.

Updated ODE function:

function zdot = odess(t,z)
% defining global variables
global r_p r_g I_p I_g m_p m_g k_mb
%%state space equations are defined below:
zdot(2) = ((r_p*k_mb)/I_p)*(-r_p*z(1)+r_g*z(3)+z(5)-z(7));  % Eq 1ss
zdot(6) = (k_mb/m_p)*(r_p*z(1)-r_g*z(3)-z(5)+z(7));         % Eq 2ss
zdot(4) = ((r_g*k_mb)/I_g)*(-r_g*z(3)+r_p*z(1)+z(7)-z(5));  % Eq 3ss
zdot(8) = (k_mb/m_g)*(-r_p*z(1)+r_g*z(3)-z(7)+z(5));        % Eq 4ss

zdot(1) = z(2); % Eq 5ss
zdot(3) = z(4); % Eq 6ss
zdot(5) = z(6); % Eq 7ss
zdot(7) = z(8); % Eq 8ss

zdot(9) = z(1); % Angular displacement?.
% Outputs
zdot = zdot';

Updated main script:

% Mainscript for vibration analysis (16/01/2017).
clear all;
clc;
%%Defining global variables
global r_p r_g I_p I_g m_p m_g k_mb
k_mb = 5;  % torsional stiffness.
r_p = 10;   % pinion pitch(?) radius
r_g = 15;   % gear pitch(?) radius
I_p = 100;  % pinion inertia.
I_g = 200;  % gear inertia.
m_p = 50;   % pinion mass
m_g = 100;  % gear mass
%==========================================================================
%%initial coditions for ode solver
%==========================================================================
z0 = [1,0,0,0,0,0,0,0,0]; % zero initial conditions.
t_step = [0:1:10];        % timestep 0 to 5 seconds.
opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3);  % conservative tolerances.
[t,z] = ode45(@odess, t_step, z0, opts);
%==========================================================================
%%Screen Output and Plotting 
%==========================================================================
[t,z(:,:)]
p_velocity = t./z(:,9);
subplot(3,1,1)
plot(t,z(:,9));
subplot(3,1,2)
plot(t,p_velocity);

[Torsten]

If you set up your code as above, the variable z(9) is z(1), integrated over time, i.e.

z_9(t)=integral_{0}^{t} z_1(t) dt.

If you want to output z(1), just use your previous code without any changes and use the command

   plot (t,z(:,1));

after the call to ode45, as I already suggested.

Your assertion

z(1) isn't included in the output matrix "z" from [t,z] = ode45(@odess,t_step,z0, opts), i think what is included are only zdots.

is simply wrong.

Best wishes

Torsten.

[Jan]

@Tyler: An integrator as ODE45 gets the start position (or "state vector") as input, than calls the function to be integrated, which replies the derivatives and accumulates them. For each successful time step, the new position is replied. Therefore the output of the integrator contains the positions (state varaibales) and the first element is its first component.

You can simply check this manually: Computer your odess(t,z) for the inital time and start vector and check the first row of the output of the integrator.

Answer 2

[Walter Roberson]

No, that cannot be done. See though

https://www.mathworks.com/matlabcentral/answers/311625-ode45-how-to-print-a-variable#answer_242828