fzero error in fa value

2 views (last 30 days)
Matias Calzada Gianero
Matias Calzada Gianero on 6 Dec 2020
Answered: Alan Weiss on 7 Dec 2020
Hi Guys!
I´m trying to implement a code that is a function which 3 parameters as arguments. It´s a model for Comparmental model study. I get the following error:
Error in fzero (line 226)
fa = localFirstFcnEval(FunFcn,FunFcnIn,a,varargin{:});
And as fzero fails, that leads to the next error inherent to the function
Error in cortinez (line 46)
ke0 = fzero(@getKe0, ike0, options);
The only arguments that I had to put on the function are age, weight and a exponent. I will let the code for anyone that maybe can find the error. But It´s an internal error, so I´m a bit lost. Thanks in advance!
function [sys] = cortinez(Age,Weight,stepSize)
% Implements the model of Cortinez using table S2a for the pooled data
% Cortínez LI, Anderson BJ, Penna A, Olivares L, Muñoz HR, Holford NH, Struys MM, Sepulveda P.
% Influence of obesity on propofol pharmacokinetics: derivation of a pharmacokinetic model.
% Br J Anaesth. 2010;105:448-56.
% Standard volumes (Liters for 70 kg patient)
V1std = 4.48;
V2std = 21.2;
V3std = 237;
%Standard clearances (Liters/min for 70 kg patient)
CL1std = 1.92;
Q2std = 1.45;
Q3std = 0.86;
%Age correction factors
SLV2 = -0.0164;
SLQ2 = -0.0153;
wa = Weight/70;
wap = wa^0.75;
aa = Age - 50;
V1 = V1std * wa;
V2 = V2std * wa * exp(SLV2*aa);
V3 = V3std * wa;
volume = [V1 V2 V3];
CL1 = CL1std* wap;
CL2 = Q2std* wap * exp(SLQ2*aa);
CL3 = Q3std*wap;
clearance = [CL1 CL2 CL3];
% Solve for the value of ke0 which yields a time to peak effect of 1.4
% minutes. We search in the range of .1 to 2.0
tpeak = 1.6;
ike0 = [.1;2];
t=0:.0005:3;
u=zeros(size(t));
u(1)=1;
options = optimset('Display','off'); % Turn off Display
ke0 = fzero(@getKe0, ike0, options);
%disp(ke0);
function [error] = getKe0(tke0)
[sys] = mam2ss(clearance,volume,tke0);
[y,t] = lsim(sys,u,t);
s = lsiminfo(y,t,0); % final value is 0
error = tpeak - s.MaxTime;
end
sys = mam2ss(clearance,volume, ke0);
if (nargin==3)
sys=c2d(sys, stepSize,'zoh');
end
end
function [sys] = mam2ss(clearance,volume, ke0)
% returns a state space model from vectors of clearance
% and compartment volumes. All states are drug amount
% in the compartment; to get concentration, divide by compartment
% volume in the C matrix. Note that this code can handle an
% arbitrary number of compartments, but assumes the first compartment is
% the central compartment (mam is short for mamillary).
% The effect site compartment is represented as a compartment of negligible
% volume adjoined to the central compartment.
n = length(clearance);
if (nargin == 2)
k1 = clearance./volume(1);
k2 = clearance(2:n)./volume(2:n);
A = [-sum(k1) k2;transpose(k1(2:n)) -diag(k2)];
B = [1;zeros(n-1,1)]; % Direct addition of drug into the central compartment
C = [1/volume(1) zeros(1,n-1)]; % Observation of the central compartment concentration.
else
EFFECT_VOL_FACTOR=10000; % ratio of central compartment vol to effect
k1 = [clearance./volume(1) ke0/EFFECT_VOL_FACTOR];
k2 = [clearance(2:n)./volume(2:n) ke0];
A = [-sum(k1) k2;transpose(k1(2:n+1)) -diag(k2)];
B = [1;zeros(n,1)]; % Direct addition of drug into the central compartment
C = [zeros(1,n) EFFECT_VOL_FACTOR/volume(1)]; % Observation of the effect site compartment concentration.
end
C = C./1000; % All compartment volumes are in L, convert to ml
D=0;
sys = ss(A,B,C,D);
end

Answers (1)

Alan Weiss
Alan Weiss on 7 Dec 2020
Next time, please give the full error thrown by the function. As it is, I don't know what the error is.
I guess that the error is that your function @getKe0 does not differ in sign at the two endpoints ike0. But that is just a guess because I don't know what the error is.
Alan Weiss
MATLAB mathematical toolbox documentation

Categories

Find more on Statics and Dynamics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!