integration of a multiple anonymous function

7 views (last 30 days)
NICOLA on 17 Oct 2014
Commented: Mike Hosea on 20 Oct 2014
Dear All,
I have a function which depend on a parameter of an integral, like the following one
g(c) = integral(x^2 + c*x + 1)
where the integration range is [0,1] and the integration variable is x. In matlab this function can be defined as a multiple anonymous function
g = @(c) (integral(@(x) (x.^2 + c*x + 1),0,1));
Now, what I need to do is to integrate a function of g(c), let's say the fourth power, over c in the same range [0,1]. Hence I wrote the code
I = integral(@(c) g(c).^4,0,1)
but it doesn't work and the reason seems to be the inner product 'c*x'. Indeed I have got the same error even if I simply do the integral of g: I = integral(g,0,1) I can simply sample the c axis and use trapz routine instead of integral or implement other quadrature schemes, but still I would like to figure out why it does not work in this way and if I can overcome the problem. Any suggestion?
Thanks Nicola

Mike Hosea on 18 Oct 2014
Edited: Mike Hosea on 18 Oct 2014
So if your function g works with a scalar value of c, then you need to vectorize it. You'll have the same problem if you do something like
x = linspace(0,1);
plot(x,g(x))
The easiest way to vectorize g is with arrayfun:
gv = @(c)arrayfun(g,c);
then something like @(c)gv(c).^4 will be vectorized properly for plotting, integrating, or whatever.

NICOLA on 20 Oct 2014
Thank you, it works perfectly. Unfortunately, this does not work (or takes very long time) when the function g depends on more parameters and the integral through which g is defined is two dimensional integral. Specifically the integral that I'm trying to solve numerically has the following shape
I = \int dx dy dz p(x,y,z) [g(x,y,z)]^4
where p(x,y,z) is a well defined density function and g(x,y,z) = \int dx' dy' f(x',y';x,y,z)
Now g is a function of three variables x,y and z and this function is defined by means of a double integration. If I apply the same procedure as above I should write
g = @(x,y,z) integral2(@(x',y') f(x',y';x,y,z),-lim,lim,-lim,lim);
gv = @(x,y,z) arrayfun(g,x,y,z);
I = integral3(@(x,y,z) p(x,y,z).*gv(x,y,z).^4,-lim,lim,-lim,lim,-lim,lim);
where the integrations are all performed over a box with linear dimension 2*lim .
Is the code right? Do you have any tips to speed up computation? Thanks in advance.
Mike Hosea on 20 Oct 2014
For technical reasons, nested integration of smooth functions often results in extra accuracy, but nothing comes with out cost. You might want to loosen the tolerances a bit, say add
'AbsTol',1e-5,'RelTol',1e-3
to the integral2 and integral3 calls. Of course, if you are nesting integral3 and integral2, it means that you are ultimately doing the equivalent work of a 5-D integral. Nested adaptive quadrature is often the first-attempted method, since it is easy to try, but you might need to resort to sparse grid or Monte Carlo methods for better speed.