Definite integration of a function

Hi All,

I am having a set of functions of time t. Some of the newer functions are defined in terms of the earlier functions. First, there are three exponential functions h(t), z(t) and r(t) that are defined as piece wise. Then there are muvI(t), Nv_adjusted(n), del_Nv_adjusted(t) and del_landings_averted(t), that are defined using the earlier functions. Then finally I want to find the definite integral of del_landings_averted to get the the value of "landings_averted". The limits of integration are "gamma" and "gamma+epsilon". When I use the quad function it gives me the error:

>> landings_averted=quad('del_landings_averted(1)',gamma,gamma+epsilon);
??? Error using ==> inlineeval at 15
Error in inline expression ==> del_landings_averted(1)
 Undefined function or method 'del_landings_averted' for input arguments
 of type 'double'.

Error in ==> inline.subsref at 27
    INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr,
    INLINE_OBJ_.expr);

Error in ==> quad at 77
y = f(x, varargin{:});

I am not able to understand how to get the value of the definite integral. I am giving the code below:

% generating the values of exponential decay of lethal effect indoors

tfinal=365; % user gives the number of days for which the simulation will be run
epsilon=90; % user gives the number of days for which the insecticide effect remains

gamma=32; % spraying is done of 1st feb of the year
Ct0=0.054;b1=-0.0260;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        h(t)= Ct0*exp(b1*tau);
    else
        h(t)=0;
    end
end

% generating the values of exponential decay of lethal effect outdoors
Ct0=0.054;b2=-0.0610;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        z(t)= Ct0*exp(b2*tau);
    else
        z(t)=0;
    end
end

% generating the values of exponential decay of repellent effect 
Ct0r=0.63;b3=-0.0130;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        r(t)= Ct0*exp(b3*tau);
    else
        r(t)=0;
    end
end

% net vector death rate with intervention
muv=0.071;
percent_landing=0.65;
Q=0.2554;
H=7933615;
C=5392890;
Hs=3000000;
Cs=2000000;

for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
        tau=t-gamma;
        muvI(t)= Q*(Hs/H)*(1-percent_landing)*r(t)*muv + Q*(Hs/H)*(1-percent_landing)*(1-r(t))*(h(t)+muv)+ Q*(Hs/H)*(percent_landing)*(h(t)+muv)+ Q*(1-Hs/H)*muv + (1-Q)*(Cs/C)*(1-percent_landing)*r(t)*muv + (1-Q)*(Cs/C)*(1-percent_landing)*(1-r(t))*(z(t)+muv)+(1-Q)* (Cs/C)*(percent_landing)* (z(t)+muv)+(1-Q)*(1-Cs/C)*muv;
    else 
        muvI(t)= muv;
    end
end 

% Adjusted sand fly population gamma days after insecticide application
% Note: the values are generated for discrete days of the years, denoted by n 
Nv=10000; % assumed the vector population as constant before spray
for n=1:tfinal
    if(n>=1 && n<=(gamma-1))
        Nv_adjusted(n)=Nv;
    else
        Nv_adjusted(n)=Nv_adjusted(n-1)*(1-muvI(n)+muv);
    end 
end

plot(h)
plot(z)
plot(r)
plot(muvI)
plot(Nv_adjusted)

% rate of change of adjusted vector population
for t=1:tfinal
    if (t>=gamma && t<=tfinal)
        del_Nv_adjusted(t)=(muvI(t)-muv)*Nv_adjusted(t);
    else
        del_Nv_adjusted(t)=0;
    end
end

% rate of change of landings averted
h_bar=5.4;
c_bar=4;

for t=1:tfinal
    if (t>=gamma && t<=tfinal)
        del_landings_averted(t)=Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+ (1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar);
    else
        del_landings_averted(t)=0;
    end
end

plot(del_Nv_adjusted)
plot(del_landings_averted)

% Need to write code for landings_averted, L(t) using integration

% len = quad(@hcurve,gamma,gamma+epsilon);
del_landings_averted(t)=Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+ (1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar);
landings_averted=quad('Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+ (1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar)',gamma,gamma+epsilon);

Please advise on how to resolve this issue.

Thanks in advance.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 1

[Walter Roberson]

Your error message is for a line that is not shown in the code you posted.

quad() is documented to only accept function handles for the function. You are attempting to pass a string instead. There is no documented behavior for such a thing.

It appears that quad() is attempting to convert the string you give in to an inline function. When it does that, it does not substitute in the current value of any variable. The inline function would be evaluated in the base workspace.

When your string 'del_landings_averted(1)' was examined in the base workspace, it saw there was no variable in the base workspace named 'del_landings_averted' and so it attempted to call a function named 'del_landings_averted' with the argument 1. No such function name could be found from the base workspace, so it generated the error you saw.

If we imagine that 'del_landings_averted(1)' had been evaluated to the value of del_landings_averted(1) at the time that quad() was invoked, then the result would have been a constant rather than a formula of some sort. Did you want to integrate the constant over a definite integral range? If so you know the result would be 0.

Note that you are not defining a function named del_landings_averted in your code: you are creating a numeric array defined over the indices 1 to tfinal. If you want to do numeric integration of that array, you could use trapz() but not quad().

Your quad() code for landings_averted should probably be replaced by

landings_averted = quad(@(t) Q .* del_Nv_adjusted(t) .* (Hs./H) .* (1-percent_landing) .* (r(t)) .* (t-gamma) ./ (Hs .* h_bar) + (1-Q) .* del_Nv_adjusted(t) .* (Cs./C) .* (1-percent_landing) .* (r(t)) .* (t-gamma) ./ (Cs .* c_bar), gamma, gamma+epsilon);

I wrote this to take in to account that quad() will pass a vector in to the function.

You still need to do more work on this, though, as the values passed in will seldom be integral, and that is going to mess up del_Nv_adjusted(t) which is access the vector del_Nv_adjusted at index t; likewise r(t) would access vector r at index t. Which would fail for t not an integral value. You are probably going to have to write del_Nv_adjusted and r and del_landings_averted as functions (or anonymous functions.)

I suggest you read http://www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html and http://www.mathworks.com/help/techdoc/math/bsgprpq-5.html

[Kaushik]

Hi Walter,

Thanks for your inputs.

First I executed the t=1:365; and

then when I execute the command that you have given above, it gives the below error:

??? Subscript indices must either be real positive integers or logicals.

Error in ==>

@(t)Q.*del_Nv_adjusted(t).*(Hs./H).*(1-percent_landing).*(r(t)).*(t-gamma)./(Hs.*h_bar)+(1-Q).*del_Nv_adjusted(t).*(Cs./C).*(1-percent_landing).*(r(t)).*(t-gamma)./(Cs.*c_bar)

Error in ==> quad at 77

y = f(x, varargin{:});

All variables here are positive integers. I do not get why I am getting this error.

Next I wrote a separate function called xy

function f = xy(t)

global Q gamma percent_landing H C Hs Cs h_bar c_bar

f=Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+ (1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar);

and then called the quad function again using this command in the script.

landings_averted = quad(@xy, gamma, gamma+epsilon);

It now gives the error,

??? Undefined function or method 'del_Nv_adjusted' for input arguments of type 'double'.

Error in ==> xy at 5

f=Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+

(1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar);

Error in ==> quad at 77

y = f(x, varargin{:});

Can you please guide me.

Thanks in advance,

Kaushik.

[Walter Roberson]

As I wrote near the bottom:

-----

You still need to do more work on this, though, as the values passed in will seldom be integral, and that is going to mess up del_Nv_adjusted(t) which is access the vector del_Nv_adjusted at index t; likewise r(t) would access vector r at index t. Which would fail for t not an integral value. You are probably going to have to write del_Nv_adjusted and r and del_landings_averted as functions (or anonymous functions.)

-----

And indeed, your code failed because it attempts to access r(t) when t is not integral. The cure, as I noted, is to write del_Nv_adjusted and r and del_landings_averted as functions (or anonymous functions.)