Get variable out of ode 45

Hi,

I have a function that goes into ode45. In this function multiple variables are calculated to finaly calculate the variable x which is the needed one. This works oke, but to review I would also like to see the other variables calculated, like q12, q23 and q34. I would not only like to see this variables. I also like to make some plots of those values against time, so disp(...) is not a good option. Please if you know how to extract those variables from the function. Let me know.

The function is used to compute water content in different soil layers. This is based on the flow q. Therefore q needs to be calculated, but I would like to know q. Here is the function given where the ode45 is used:

% Calculate content per layer
fh = @(time,x)fflow(time, x, tu, p, Ks, theta_r, theta_s, lambda, n, alpha);    % Create function handle to put into ODE-fucntion
[time, xt] = ode45(fh, time, x0);                                            % Compute water content 'x' per layer at time t

The different q values are computed in fflow, part of the code is given for the computation of q12 but for other q values it is the same:

if (TF(1) == 1) && (TF(2) == 0)  || ((TF(1) == 1) && (TF(2) == 1))                 % First layer saturated, second layer not: 
    q12 = maxFlowq(Ks(1),lambda(1), n(1), 30);  % Maximum flow rate for soil type
elseif ((TF(1) == 0) && (TF(2) == 1))          % First layer not saturated, second layer saturated
    q12 = 0;
else                                            % First and second layer not saturated
    q12 =   -K(1)*((-abs(diff(1))/p(1))+1);
end

% Compute change in water content theta(x)
dx1dt = (qin-q12)/p(1);                % Water content difference layer 1 
dx2dt = (q12/p(1)-q23/p(2));           % Water content difference layer 2 
dx3dt = (q23/p(2)-q34/p(3));           % Water content difference layer 3

dxdt = [dx1dt; dx2dt; dx3dt];

So if possible I would like to have the values for q for every t from this function.I hope you can help me:) Kind regards, Jits

Answer 1

[michio]

One way is to use OutputFcn. The following entry could be of use.

How do I extract extra parameters from the ode45 ordinary differential equation solver function?

For details on this option, please see "OutputFcn" section of odeset.

[Jits Riepma]

Hi,

Thank you for the answer! I tried it. The implementation of the odeset with OutputFcn in the main code is now correct I think? :

% Calculate content per layer
fh = @(time,x)fflow(time, x, tu, p, Ks, theta_r, theta_s, lambda, n, alpha);    % Create function handle to put into ODE-fucntion
options = odeset('OutputFcn',@(time,x,flag) myOutputFcn(time, x, flag, tu,  p, Ks, theta_r, theta_s, lambda, n, alpha)); 
[time, xt] = ode45(fh, time, x0, options);                                            % Compute water content 'x' per layer at time t

But I have no idea how to fill in the function of "myOutputFcn". Should I fill in the three different if statements for every q at the init and the ''? I understand that at the 'done' the values of computed are assigned to different variables for in the workspace. But I don't understand what happends before .

[michio]

After initialization, q1234 is computed at each time step and added to q_save, ie. the size of q_save increases as time steps and ends up to be Mx3, where M is the number of time steps taken to solve the ODE.

An example myOutputFcn for you is the following. Hope it resolves your issues.

function status = myOutputFcn(time, x, flag, tu,  p, Ks, theta_r, theta_s, lambda, n, alpha)
% OutputFcn sample
%
% Calculate q's (q12, q23, q34) here
% Suppose these are scalar values and let's store into one variable q1234
% ex: q1234 = [q12, q23, q34];
% 
%
persistent q_save
switch flag
    case 'init'
        q_save = q1234;
    case ''
        q_save = [q_save; q1234];
    case 'done' % when it's done
        assignin('base','q_save',q_save); % get the data to the workspace.
end
status = 0;

[Jits Riepma]

Hi,

Thank you! This kind of works! The only thing is that it only saves 13 values for q. This because it takes time steps bigger than 1, if I check the proces with the debugger the values for t some times are like: [3 4 5 6 7]. Then it only saves the q1234 for t = 7 Is there some way to decrease the time step t, so that it takes only steps of 1?

[michio]

I see... I would suggest using tspan (interval of integration) with a longer vector of the form [t0,t1,t2,...,tf] instead of two elements. myOutoutFcn will be evaluated at one time step at a time.

In your code, execute

[time, xt] = ode45(fh, tspan, x0, options);

with tspan with a vector,

 tspan = linspace(t0,tf,100);

for example. As stated in the documentation page of ode45

"If tspan has two elements, [t0 tf], then the solver returns the solution evaluated at each internal integration step within the interval. If tspan contains more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points."

Answer 2

[Torsten]

Try the suggestion under

https://de.mathworks.com/matlabcentral/answers/92701-how-do-i-pass-out-extra-parameters-using-ode23-or-ode45-from-the-matlab-ode-Suite

Best wishes

Torsten.

Answer 3

[Jan]

Do not integrate dicontinous functions with Matlab's ODE integrators. They are not designed to handle this correctly and the result can be wrong, dominated by rounding errors or if you are lucky the integration stops with an error:

http://www.mathworks.com/matlabcentral/answers/59582-ode45-errer-input-argument-y-is-undefined#answer_72047