Evaluate numerical Integral by generating function

I want to evaluate numerically a double integral of a function.

The function is generated from a series of steps. My doubt is how to handle the variable, while generating the function.

The function is generated from a series of many steps. For example, I am giving a small version of it, where I tried by using handles but I am getting errors. Kindly let me know the correct way of doing it.

ss = 10^-8;
sd = 10^-11;

ber_link = 0;
tn = 6;
ma = [1,1;1,-1;1,0;-1,1;-1,-1;-1,0];
p = [.9,.01,.01,.01,.01,.01];

for j = 1:tn
        x = 0;
        x = x+ma(j,2)*ss;
        if ma(j,1) == 1
            a = @(ln1, ln2) (ss*ln1+x*ln2)/sd;
        else
            a = @(ln1, ln2) (ss*ln1-x*ln2)/sd;
        end

          cpd = @(ln1, ln2) qfunc(a);
          pro = @(ln1, ln2) cpd*p(j);

          bl = @(ln1, ln2) bl+pro;

end    

bl = @(ln1,ln2) bl.*exp((-(ln1-3).^2)/2*36).*exp((-(ln1-3).^2)/2*36);
int = integral2(bl, 1, Inf, 1, Inf)

Answer 1

[Mike Hosea]

This

bl = @(ln1, ln2) bl+pro;

is syntactically wrong. The operation + is not defined on function handles. You need to include the arguments. The statement should read

bl = @(ln1,ln2) bl(ln1,ln2) + pro(ln1,ln2);

You'll need to make similar changes in the definition of cpd and pro. Try some simple examples before jumping to the real thing to make sure you understand how the language works first. Something like:

>> f = @(x)ones(size(x));
>> f = @(x)f(x).*x + 1; % now f(x) = x + 1
>> f = @(x)f(x).*x + 1; % now f(x) = x^2 + x + 1
>> f = @(x)f(x).*x + 1; % now f(x) = x^3 + x^2 + x + 1
>> f(1)
ans =
     4

The next problem you will encounter here is that you are using matrix multiplication * and matrix division (linear system solves) / when what you really mean is element-wise multiplication .* and element-wise division ./. This will become important when you call integral2 because integral2 requires that your function be able to accept array input and operate element-wise. You are doing this at the last step of defining bl, but it is just as important in the expressions in the loop as when finishing off the definition. Just FYI, if you have a function that can't take an array, you can always integrate

   @(x,y)arrayfun(f,x,y)

instead of integrating just f. ARRAYFUN takes care of calling f element-wise with only scalar inputs. I don't think you'll need that here, however, as it looks like everything can be done element-wise.

[Stephen23]

Super answer, very clearly explained. I would accept this if I was the question asker!

[Guillaume]

There also seem to be some confusion between function handles and variables. This line:

bl = @(ln1,ln2) bl.*exp((-(ln1-3).^2)/2*36).*exp((-(ln1-3).^2)/2*36);

declares a function handle bl and use a variable bl within the expression.

Also:

bl = @(ln1,ln2) bl(ln1,ln2) + pro(ln1,ln2);

Is never going to work properly. There's a bl too many in this expression.

[Sri]

Thank you for the reply.

@Guillaume, I wanted to use

bl = @(ln1,ln2) bl(ln1,ln2) + pro(ln1,ln2);

as I wanted to add iteratively to get final function for integration. As you pointed out, it gave me error and I used sum and modified the code as given below.

ss = 10^-8;
sd = 10^-11;

ber_link = 0;
tn = 6;
ma = [1,1;1,-1;1,0;-1,1;-1,-1;-1,0];
p = [.9,.01,.01,.01,.01,.01];

for j = 1:tn
        x = 0;
        x = x+ma(j,2)*ss;
        if ma(j,1) == 1
            a = @(ln1, ln2) (ss*ln1+x*ln2)/sd;
        else
            a = @(ln1, ln2) (ss*ln1-x*ln2)/sd;
        end

          cpd = @(ln1, ln2) qfunc(a(ln1,ln2));
          pro = @(ln1, ln2) cpd(ln1, ln2)*p(j);

          bl = @(ln1,ln2) sum(pro(ln1, ln2));

end    

func= @(ln1,ln2) bl(ln1, ln2).*exp((-(ln1-3).^2)/2*36).*exp((-(ln1-3).^2)/2*36);
int1= integral2(func,1,Inf,1,Inf);

Could someone let me know, if it is correct?

I tried looking into arrayfun, but it is not clear. Could you kindly elaborate little more on it.

[Sri]

Have another doubt. In the above code, if the following code changes as given below

if ma(j,2) == 0
 a = @(ln1) (ss*ln1)/sd;
 else
 a = @(ln1, ln2) (ss*ln1-x*ln2)/sd;
 end
cpd = @(ln1, ln2) qfunc(a(ln1,ln2));
pro = @(ln1, ln2) a(ln1, ln2)*p(j);
bl = @(ln1,ln2) sum(pro(ln1, ln2));

How to handle 'cpd', 'a' and 'pro' in this case, as in one case they are function of ln1 and in other case, they are function of two variables (ln1, ln2).

I tried in the following way

if ma(j,2) == 0
a1 = @(ln1) qfunc((ss*ln1)/sd)*p(j);
else
a2 = @(ln1, ln2) qfunc((ss*ln1-x*ln2)/sd)*p(j);
end

bl = @(ln1,ln2) sum(a1(ln1)+a2(ln1, ln2));

Is this fine to do this way?

[Mike Hosea]

I think I would just make it a function of two arguments always. You're really just defining the maximum number of arguments that the function will accept, but you don't have to reference the second argument.

>> f = @(x,y)cos(x)
f = 
    @(x,y)cos(x)
>> f(pi/4)
ans =
   0.707106781186548
>> f(pi/4,'How is the weather over there?')
ans =
   0.707106781186548
>> f(pi/4,2,3)
Error using @(x,y)cos(x)
Too many input arguments.

BTW, I think Guillaume did not quite understand what you were trying to do. I have to say, this is one of the most advanced uses of function handles, taking advantage of the way they "snapshot" the referenced variables at the time of their creation, that I have come across, but this is the sort of thing that is a possibility in an interpreted language. Many years ago I was a LISP programmer, and you are bringing back some warm memories. :)