How do I extract a function calculated and used inside my ode45 function?

Hello,

I'm using ode45 to produce solutions to ODE's and It works. However, I have a function inside of the code that I want to do a graph in the main.m. The system is modelled with the following equations:

Well, the graph I want to do is and versus time. Do you have any ideas how to get them?

Thank you very much in advance

[Jan]

This is easy, if you post the code of the function to be integrated: It does contain the calculation of h, so simple return it as 2nd output. Therefore the function to be integrated must accept the input as vector.

[Ignacio López]

Hello Jan, thanks for your answer. I didn't understand very well what you propose...

I tried to add 2 aditional equations, using this idea:

Likes this, Ode45 would retur and versus time., but from a numerical point of view is a bit tricky, because when I tried to do the derivative of these functions I used the function "diff" divided by a small number, but this expression does not have the same length of the other vectors...

[Jan]

I cannot follow you. Please post the current code. As far as I understand, you calculate h and y already inside the function to be integrated. Then it is very easy to export them. Describing the method would take a lot of time, but implementing it is most like just some key strokes, when the function is available.So please be so kind and show us the code.

[Jan]

[MOVED from section for answers] Ignacio López wrote:

Hello Jan, once again thanks for answering :)

So the first equation the biomass balance and the second one the oxygen generation (I cannot post the whole code, because some parameters "p" are in the main.m)

dy(1) = u_Net*y(1)-Ru*(u/(Ag.^2))*y(1); 
dy(2) = k_O.*mu_bar(I,y(1),p)*y(1) - (kL_O2*(y(2) - O2_star) - ...
        RO2.*(u/(Ag.^2))*y(1)); 

The curve of oxygen reaches a maximal value, and I would like to show that at first the generation term (the one of the right in equation 2) is stronger and consumption, and after reaching the maximum, the behaviour is on the other way around; consumption stronger than generation. And indeed when I do the graph dO2/dt versus time, I obtained the expected results, but I wanto show that is due to the reason I explained before.

Do you understand me? Maybe I'm not being clear...

Thanks again!

[Jan]

I do not understand the equation or the description. But understanding the meaning is not needed to solve the problem. You've asked for getting the values of h and y. I do not see, where the values of thses functions are calculated. Therefore it is still not clear to me, what you are asking for.

As far as I understand, you need something like this:

[t, y] = ode45(@(t, y) fcn(t, y, p), [0, 1], [0, 0]);

function dy = fcn(t, y, p)
a  = t * p(1);  % or what ever
b  = y(1) * t^2;
dy = [a * y(2); b * y(1)];
end

Now you want to obtain the values of a and b for the evalulated time steps. Then:

[t, y] = ode45(@(t, y) fcn(t, y, p), [0, 1], [0, 0]);
[dy, a, b] = fcn(t, y.', p);  % <== here a and b are requested

function [dy, a, b] = fcn(t, y, p)
a  = t * p(1);  % or what ever
b  = y(1, :) * t .^ 2;
dy = [a .* y(2, :); b .* y(1, :)];
end

I'm not sure about the transpositions of y. Please debug this. The point is, that the function must accept a vector as input and replies the wanted arrays as 2nd or 3rd output argument, where they are ignore during the integration, but it is easy to calculate them afterwards, when the t and y values are available.

Does this match your question?

Answer 1

[Ignacio López]

Thanks Jan for your answer!

But I can realise that maybe you are making the same mistake I was doing it. I mean, in this function

function [dy, a, b] = fcn(t, y, p)
a  = t * p(1);  % or what ever
b  = y(1, :) * t .^ 2;
dy = [a .* y(2, :); b .* y(1, :)];
end 

you are taking the integrated valures of "y" to create "a". I used your way to do it, and I obtained the same graph I did before. What I actually need are the values of [a * y(2); b * y(1)] while ode45 is running... Because I want get the behaviour of the terms which modify "dy".

Thank you very much for your time Jane!

[Jan]

@Ignacio: If you need dy, you can use the 1st output of the function directly. Remember that "while ode45 is running" means even the rejected steps, which are not included in the replied trajectory. This is rarely useful.