integration of a multiple anonymous function

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

Answer 1

[Mike Hosea]

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.

Answer 2

[NICOLA]

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]

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.