MATLAB: Get variable out of ode 45

function variable from ode45MATLABode45outputfcn

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

Best Answer

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.

Related Solutions

MATLAB: Different input per time step ODE45

function main
p = [25;25;25];           %(constant) depth of soil layers [cm]
% If time is not imported
time = (0:1:150)';     % time vector
x0 = [0.058;0.042;0.046];       % vector for initial water content values   
% Calculate water flow
fh = @(t,x)fflow(t, x, time, u, p);
[T,X] = ode45(fh, time, x0);
function dxdt = fflow(t,x,time,u,p)
qin = interp1(time,u,t);                  % Volume flux density as input (=irr-ET)
q12 = -K(1)*((-abs(diff1)/p(1))+1);       % Volume flux density layer 1 to 2 [cm3 cm-2 h-1]
q23 = -K(2)*((-abs(diff2)/p(2))+1);       % Volume flux density layer 2 to 3 [cm3 cm-2 h-1]
q34 = -K(3)*((-abs(diff3)/p(3))+1);       % Volume flux density layer 3 to 4 [cm3 cm-2 h-1]
% Compute change in water content theta (x)
dx1dt = (qin-q12)/p(1);       
dx2dt = (q12/p(1)-q23/p(2));
dx3dt = (q23/p(2)-q34/p(3));
dxdt = [dx1dt; dx2dt; dx3dt];

Best wishes

Torsten.

MATLAB: How to save a variable computed inside ODE

assignin is a bad programming technique. Producing variables remotely in another workspace is not clean. A better method is to let the integration run regularily and request the wanted values afterwards:

[t, y] = ode45(@fcn, tSpan, y0);
function dy = fcn(t, y)
p  = sin(y(2));   % <- The wanted parameter

dy = [2 * t; p / y(1)];
end

If this is the function to be integrated, it must be modified such, that it accepts vectors as input and replies p as output:

function [dy, p] = fcn(t, y)
p  = sin(y(:, 2));   % <- The wanted parameter
dy = [2 .* t; p ./ y(:, 1)];
end

Then:

[t, y] = ode45(@fcn, tSpan, y0);
[~, p] = fcn(t.', y.');

This is much easier than letting the output function create variables remotely.

For some cases the output function is the best choice. Then provide the value as output instead of assignin

function status = myOutputFcn(time, x, flag)
persistent q_save
status = 0;
switch flag
    case 'init'
        q_save = q;
    case ''
        q_save = [q_save; q];
    case 'done' % when it's done
        status = q_save;
        q_save = [];
end
end

But now to the actual problem: Unfortunately the ODE integrators reply different time steps, if the output is obtained as a struct. You did not show, how you call the integrator, so the most important detail is missing. (In opposite to this, the code or the output function is not a part of the problem.)

If you want to call the integrator in the form sol = ode45(___), and use a tspan defined as vector also, you have to modify ode45 such that the Output function is called regularily as for the [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) form. I consider this a bug and sent a bug report to MathWorks some years ago.

You wrote:

"I have also tried declaring tspan in two ways"

  tspan = linspace(t0,tf,100)  
  tspan = [t0 tf]

And what did you observe as effect? How did you call which integrator?

Related Question