BYOM function call_deri.m (calculates the model output)
Syntax: [Xout,TE,Xout2,zvd] = call_deri(t,par,X0v,glo)
This function calls the ODE solver to solve the system of differential equations specified in derivatives.m, or the explicit function(s) in simplefun.m. As input, it gets:
- t : the time vector
- par : the parameter structure
- X0v : a vector with initial states and one concentration (scenario number)
- glo : the structure with various types of information
The output Xout provides a matrix with time in rows, and states in columns. This function calls derivatives.m. The optional output TE is the time at which an event takes place (specified using the events function, set it up to catch discontinuities). If you want to use parameters that are (or influence) initial states, they have to be included in this function.
- Author: Tjalling Jager
- Date: November 2021
- Web support: http://www.debtox.info/byom.html
- Back to getting started or byom_my_first.m
Copyright (c) 2012-2021, Tjalling Jager, all rights reserved. This source code is licensed under the MIT-style license found in the LICENSE.txt file in the root directory of BYOM.
Contents
Start
The function statement, and make the global structure glo known here (we don't use zvd in this demo).
function [Xout,TE,Xout2,zvd] = call_deri(t,par,X0v,glo)
Xout2 = []; % additional uni-variate output, not used in this case zvd = []; % additional zero-variate output, not used in this case % Note: if these options are not used, these variables must be defined as % empty as they are outputs of this function.
Initial settings
This part extracts optional settings for the ODE solver that can be set in the main script (defaults are set in prelim_checks). The useode option decides whether to calculate the model results using the ODEs in derivatives.m, or an analytical solution that has to be provided in simplefun.m. Using eventson=1 turns on the events handling. Also modify the sub-function at the bottom of this function! See the example BYOM files for more information.
useode = glo.useode; % calculate model using ODE solver (1) or analytical solution (2) eventson = glo.eventson; % events function on (1) or off (0) stiff = glo.stiff; % set to 1 to use a stiff solver instead of the standard one % Unpack the vector X0v, which is X0mat for one scenario X0 = X0v(2:end); % these are the intitial states for a scenario
Calculations
This part calls the ODE solver (or the explicit model in simplefun.m) to calculate the output (the value of the state variables over time). There is generally no need to modify this part. The solver ode45 generally works well. For stiff problems, the solver might become very slow; you can try ode15s instead.
c = X0v(1); % the concentration (or scenario number) t = t(:); % force t to be a row vector (needed when useode=0) TE = 0; % dummy for time of events if useode == 1 % use the ODE solver to calculate the solution % specify options for the ODE solver; feel free to change the tolerances, % if you know what you're doing (for some problems, it is better to set % them much tighter, e.g., both to 1e-9) reltol = 1e-4; abstol = 1e-7; options = odeset; % start with default options if eventson == 1 options = odeset(options,'Events',@eventsfun,'RelTol',reltol,'AbsTol',abstol); % add an events function and tigher tolerances else options = odeset(options,'RelTol',reltol,'AbsTol',abstol); % only specify tightened tolerances end % call the ODE solver (use ode15s for stiff problems, escpecially for pulses) if isempty(options.Events) % if no events function is specified ... switch stiff case 0 [~,Xout] = ode45(@derivatives,t,X0,options,par,c,glo); case 1 [~,Xout] = ode113(@derivatives,t,X0,options,par,c,glo); case 2 [~,Xout] = ode15s(@derivatives,t,X0,options,par,c,glo); end else % with an events functions ... additional output arguments for events: % TE catches the time of an event, YE the states at the event, and IE the number of the event switch stiff case 0 [~,Xout,TE,YE,IE] = ode45(@derivatives,t,X0,options,par,c,glo); case 1 [~,Xout,TE,YE,IE] = ode113(@derivatives,t,X0,options,par,c,glo); case 2 [~,Xout,TE,YE,IE] = ode15s(@derivatives,t,X0,options,par,c,glo); end end else % alternatively, use an explicit function provided in simplefun! Xout = simplefun(t,X0,par,c); end if isempty(TE) || all(TE == 0) % if there is no event caught TE = +inf; % return infinity end
Output mapping
The matrix Xout contains a row for each state variable. It can be mapped to the data. If you need to transform the model values to match the data, do it here.
% Xout(:,1) = Xout(:,1).^3; % e.g., do something on first state variable, like cube it ...
Events function
Modify this part of the code if eventson=1. This subfunction catches the 'events': in this case, it looks for the start of depuration at t=7. This is rather trivial, as we already know exactly when this event happens. However, it is useful as the ODE solver will now stop and restart at t=7 and will be hampered less by the sudden shift in states at this discontinuity. A better solution would be to run the solver in two runs, first for t=0-7 and then from t=7-end. However, for this example, we won't go through that trouble.
Note that the eventsfun has the same inputs, in the same sequence, as derivatives.m.
function [value,isterminal,direction] = eventsfun(t,X,par,c,glo) nevents = 1; % number of events that we try to catch value = zeros(nevents,1); % initialise with zeros value(1) = t - 7; % thing to follow is when depuration starts at t=7 isterminal = zeros(nevents,1); % do NOT stop the solver at an event direction = zeros(nevents,1); % catch ALL zero crossing when function is increasing or decreasing