ODE - Save variables for every time step

Question

sko 2021-3-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/766021-ode-save-variables-for-every-time-step

评论： Bjorn Gustavsson 2021-4-22

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bjorn Gustavsson 2021-3-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/766021-ode-save-variables-for-every-time-step#answer_641771

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Bjorn Gustavsson 2021-3-8

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

sko 2021-3-8

在 MATLAB Online 中打开

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 2021-4-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/766021-ode-save-variables-for-every-time-step#answer_681445

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Walter Roberson 2021-4-22

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 2021-4-22

编辑：sko 2021-4-22

在 MATLAB Online 中打开

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 2021-4-22

在 MATLAB Online 中打开

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 2021-4-22

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

请先登录，再进行评论。