ODE - Save variables for every time step

Hello Guys,

I am stuck with saving time dependent variables inside a model function ... I was browing some already answered questions but none seems to help.

I have the following:

[t, f_val] = ode15s(@model,tspan,yeq);

That yields my f_val as a function of time.

But inside my model function I have also yy variable (time dependent) I would like to trace

function [f, yy] = model(t,yeq)
f= ... %some expression
yy= ... %some expression

How can I extract yy as a function of time?

Cordially

Answer 1

[Bjorn Gustavsson]

This is something that is (simplest?) done by calculating the second output once after you've integrated the ODE-system. Either in one fell swoop if your ODE-function can handle arrays for input of t and f:

[f_second_time_around_or_tilde, yy] = model(t,f_val);

Or you can quickly loop over t:

for i_t = numel(t):-1:1,
    yy(i_t,:) = model(t(i_t),f_val(i_t,:));
end

Since the odeNN-functions take cunningly adaptive steps in time we cannot trivially save all the values of yy along the steps...

HTH

[sko]

Hi,

Thank you for your contribution.

However,

[~, yy]=model(t,f_val);

fives the vaule of yy for the last time only! I want to have my output of yy's for every (converged) time step.

Best wishes

[Bjorn Gustavsson]

What with the second loop-based suggestion? That explicitly calculates yy for every step in the solution you have...

[sko]

Hello Bjorn,

I re thought my code a little bit and I think I understand my problem better, however, still without solution.

In my main code i need to 1) specify my initial conditions (initial concentrations):

%% Concentration in gas phase %%
l=1:1:ent.nz+1;
i=1:1:ent.ncg; % ncg
rx.Cg(i,l)=0; 
rx.Cg(1,l)=1; %A
%% Concentration in solid phase %%
l=1:1:ent.nz+1;
k=1:1:ent.nr;
i=1:1:ent.ncg;
rx.particle.Cp(i,k,l)=0;

2) Load those concentrations in terms of single vector (for every l,k and i I create a single vector yeq)

3) Start ODE solver

tspan=[0 10];
[t, f_val] = ode15s(@model,tspan,yeq);

The challange with my model is the fact that time derivatives I need to specify depend on concentration in space (i.e. rx.dCgdt(i,l) and rx.particle.dCpdt(i,k,l) are functions of rx.Cg(i,l) and rx.particle.Cp(i,k,l) ) Thus, my

@model 

function must first deload concentrations rx.Cg(i,l) and rx.particle.Cp(i,k,l) from loaded vector yeq, and secondly specify x.dCgdt(i,l) and rx.particle.dCpdt(i,k,l). Once time derivatives calculated as demanded by ode solver, I need to load them into single vector dydeq specified as f in the function model.

Finally,

[t, f_val] = ode15s(@model,tspan,yeq);

Yields my concentrations rx.Cg(i,l) and rx.particle.Cp(i,k,l) with time.

Answer 2

[sko]

Hello Guys,

I come once again on this general topic. I really really dont get how to code function so that I obtain multiple outputs. I will reprhase my question so that it becomes clearer.

I have an ODE solver that works just fine that I call as follows:

[time,concentration] = ode15s(@(t,yeq)model(t,yeq,dr,nrp,Deff),timespan,yeq);
c=concentration';

The function that calculates time derivatives dcdt is model function. However, inside this function I also want to calculate some others parameters regarding kinetics (not only concentrations I obtain from ODE solver) such as conversion, production selectivities that all depend on calculated concentration and are thus time dependent. Therefore, I need to correct my model function so that not only it provides me calculated concentrations but also those other parameters. Although I calculate those parameters in my model function, I dont seem to get how to extract them in my main file.

My model function is already structured so that I calculate time derivatives dcdt but also selectivity etc ...

Someone please help ....

[Bjorn Gustavsson]

Write your function something like this:

[dcdt,Par1,Par2] = model(t,C_of_t,reaction_coeffs_etc)
k1 = reaction_coeffs_etc(1);
k2 = reaction_coeffs_etc(2);
dcdt = [-k1*C_of_t(1)*C_of_t(2);
        -k1*C_of_t(1)*C_of_t(2);
        k1*(k2/(1+k2))*C_of_t(1)*C_of_t(2);
        k1*(1/(1+k2))*C_of_t(1)*C_of_t(2)]; % or whatever reactions you have...
par1 = something;
par2 = something_else;
if nargout > 1
  Par1 = par1;
end
if nargout > 2
  Par2 = par2;
end

HTH

[sko]

Hello Bjorn,

Thank you very much for your reponse. I modified my model function exactly like that. Now, where are my par1 and par2 variables stocked? I dont see them in my workspace (because they are evaluated inside model function), so how should I modify the line accordingly:

[time,concentration] = ode15s(@(t,yeq)model(t,yeq,dr,nrp,Deff),timespan,yeq);
c=concentration';

Thanks!

[Bjorn Gustavsson]

You simply loop over the solution (time, concentration) a second time after the call to ode15s:

[time,concentration] = ode15s(@(t,yeq)model(t,yeq,dr,nrp,Deff),timespan,yeq);
c=concentration';
for i_t = numel(t):-1:1,
  [yy(i_t,:),PAR1(i_t),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);
end

[sko]

Here is the error I am getting:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in test_production (line 125)
  [yy(i_t,:),PAR1(i_t),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);
  

I followed exactly your advice.

However, when I call

  model(time(i_t),concentration(i_t,:),dr,nrp,Deff)
  

I am still not getting my Par1 and Par2, although Matlab properly enters function model and calculate both Par1 and Par2!

[Walter Roberson]

Therefore, I need to correct my model function so that not only it provides me calculated concentrations but also those other parameters.

NO . It is very unlikely you want to do what you are asking for.

Each step of ode45() invokes the ode function multiple times, with different parameters -- different boundary conditions and different times. ode45() in particular invokes the ode function 6 times for each step. A prediction is made using the results of the first 5 calls, and then cross-checked with a prediction made using the result of the 6th call included. If the 4th order projection agrees within tolerance with the 5th order projection, then the step is accepted and recorded, and if the 4th disagrees with the 5th too much, then the step is rejected and another attempt is made from the same starting point but with a shorter step. In areas that are changing rapidly, it is normal for several trials to be rejected.

The ode function is called for all of those locations -- the 5 intermediate locations at the best of times, and possibly a number of more times as a step fails as it encounters a steep region and has to "slow down" to accomodate.

Furthermore, the places that ode45 emits into the outputs are not necessarily anywhere that ode45 has called the ode function for; ode45 periodically makes projections and saves the output of the projections.

So, we have a situation where the ode function gets invoked many times at locations that will not be used in the final output, plus the final outputs might be at locations that the ode function was never called for.

You are asking for the intermediate values for all of the calls to be reported back to you. Which will include at locations that ode45 decided were too steep and never tried to travel to. And which will not include at locations tha ode45 reported back for but never actually evaluated at.

You are unlikely to have any use for all that output: it tells you about places you do not need to know about, and it does not tell you about places you need to know about.

On the other hand, the strategy @Bjorn Gustavsson is recommending is very effective: it gives you back information for exactly the locations that were reported back to you by ode45(), which is much more likely to be what you want to know.

[sko]

I agree that times Bjorn suggested are those I am in fact looking for. However, even though my function is definied as:

function [dcdt, Par1, Par2] = model(t,yeq,dr,nrp,Deff)
dcdt = some expression
par1 = somme expression
par2 = some expression

if nargout > 1
  Par1 = par1;
end
if nargout > 2
  Par2 = par2;
end

end

When I iterate over solutions in my main by:

[time,concentration] = ode15s(@(t,yeq)model(t,yeq,dr,nrp,Deff),timespan,yeq);
c=concentration';
 
 
for i_t = numel(time):-1:1,
[yy(i_t,:),PAR1(i_t),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);
end

I still dont get my Par1 and Par2 appearing in my main. Moreover calling model function with

model(time(i_t),concentration(i_t,:),dr,nrp,Deff);

yields me a vector 55:1 which is I believe my dcdt (only first output of my solution).

[Walter Roberson]

When you call

for i_t = numel(time):-1:1,
[yy(i_t,:),PAR1(i_t),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);
end

is it complaining about dimension mismatches in the assignments? Could you confirm that numel(t) > 1 ? After that loop, what do you see in the code as size(PAR1) and size(PAR2) ?

[sko]

I think he is complaining regarding the fact that

model(time(i_t),concentration(i_t,:),dr,nrp,Deff)

yields a single vector (output dcdt of the function model) which is consisten with the left side

[yy(i_t,:),PAR1(i_t),PAR2(i_t,:)]

However, when I call

  [yy(i_t,:),PAR1(i_t),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);

Matlab properly enters model and properly calculates dcdt Par1 and Par2, for a given times but

model(time(i_t),concentration(i_t,:),dr,nrp,Deff)

outputs only dcdt...

I underline the fact tgat my function model is defined as:

function [dcdt, Par1, Par2] = model(t,yeq,dr,nrp,Deff)
...
end

So I dont seem to get what is wrong.

NB numel(time) is indeed >1 and represents number of time values sought.

NB2 After the loop, it enters model function and calculated dcdt (1x55) par1 (1x2) and par2 (1x2) while in model function. Then it leaves model function and resulting output from

model(time(i_t),concentration(i_t,:),dr,nrp,Deff)

is only 1x55 ...

[sko]

Found it!

for i_t = numel(time):-1:1,
  [yy(i_t,:),PAR1(i_t,:),PAR2(i_t,:)] = model(time(i_t),concentration(i_t,:),dr,nrp,Deff);
end

Thank you very very much both of you! It was a mismatch in dimensions of PAR that model produced.

[Bjorn Gustavsson]

Yes, the in solution I provided I couldn't now what the dimensions of your additional outputs were expected to be, therefore I opted to provide something that you could easily adjust to match the dimensions of your outputs (one or the other parameter might even have different dimensions at different time-steps as far as I could see...).