Pere Garau Burguera
Pere Garau Burguera

Reputation: 71

Integrating a multivariable function one variable at a time

I have a function which depends on x and y, f(x,y) (returns a scalar) and I need to integrate it so that for a particular point (x0,y0), I get a function p(theta), I basically convert it to polar, change the origin and integrate over r, so the remaining function only depends on theta:

%fun = @(x,y) Some function
fun = @(x,y) fun(x+x0,y+y0);
polarfun = @(r,theta) fun(r.*cos(theta),r.*sin(theta)).*r;
p = @(theta) integral(@(r)polarfun(r,theta),0,Inf);

After I get this function, I need to integrate it again (w.r.t theta). This doesn't work (a and b are the limits of the integral)

value = integral(p,a,b);

I need to set ArrayValued = true for it to work

value = integral(p,a,b,'ArrayValued',true);

I don't understand why I need to specify it's an ArrayValued function when it's not meant to be, also the integral takes much more time, and I need it to be calculated as fast as possible. Can you explain me why it is ArrayValued, and whether there is a way to make it without the ArrayValued set to true?

Upvotes: 1

Views: 468

Answers (1)

TroyHaskin
TroyHaskin

Reputation: 8401

Why

There are two levels as to why the ArrayValued parameter is required to be set true: function requirements and why the function has those requirements.

The requirement of the integrand function handle fun per its documentation is:

For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result,y. ... If you set the 'ArrayValued' option to true, then fun must accept a scalar and return an array of fixed size.

Per your code, the function handle p receives a vector of integration points via integral if ArrayValued is not set to true. Then, polarfun receives a vector of theta and integral(@(r)polarfun(r,theta),0,Inf) will produce a scalar value assuming no dimension mismatch happens while that inner-integral executes. And since the outer-integral passed a vector and received a scalar back, an error is thrown since those sizes don't match and ArrayValued is false by default.

Now, why are the requirements the way they are? Cleve Moler recently wrote a blog post discussing integral and its predecessors. The long-and-short is that for scalar integrands, vectorized function evaluations yield the best performance in Matlab; thus, it is used by default.

If the Mathworks wanted to, they could probably write a more sophisticated mechanism for carrying out the evaluation with vector-valued integrands through the use of bsxfun or a similar utility. However, the added complexity and subsequent required specific input requirements would more than likely be too burdensome for potentially little, if any, performance increase. So they assume an inherently scalar function that is Matlab vectorized with an option to fallback to scalar node evaluation. Additionally, your use-case showcases a huge problem since polarfun possesses a fixed vector of theta that will more than likely generate a dimension mismatch or give an outright wrong result if integral needs to refine its r integration at all.

Recommendation

With all of that being written, I think it's best to just use integral2 and let Matlab handle everything:

fun      = @(x,y) fun(x+x0,y+y0);
polarfun = @(r,theta) fun(r.*cos(theta),r.*sin(theta)).*r;
value    = integral(polarfun,0,Inf,a,b);

Upvotes: 1

Related Questions