Evaluate numerical Integral by generating function

4 views (last 30 days)
Sri on 20 Sep 2014
Edited: Mike Hosea on 21 Sep 2014
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;
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)

Mike Hosea on 20 Sep 2014
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.
Mike Hosea on 21 Sep 2014
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. :)