quadl with vector bounds

I want to do the following:

define a function like

fun_int=@(x) quadl(some_function_handle,0,x,TOL);

and then evaluate the integral at different bounds x by just calling the function in an algorithm.

However when i call fun_int with a vector argument it does not work, since quadl can't have vector bounds.

My workaround is to just use a loop, but that is slow and not very elegant. Any ideas how to achieve what I want to do without using loops?

Answer 1

[Mike Hosea]

One option is to use ARRAYFUN

>> arrayfun(@(b)integral(@(x)x.*sin(x),0,b,'abstol',1e-5,'reltol',1e-5),0:3)
ans =
         0    0.3012    1.7416    3.1111

If you don't have R2012a, use QUADGK. Though they offer faster results at lower accuracy in some cases (like my example above), QUAD and QUADL are obsolete. I think the most efficient approach might be to sort the vector b, integrate each subinterval, use cumsum to compute the individual integrals, and then unsort it.

function q = vint(f,a,b,atol,rtol)
[bsorted,idx] = sort(b);
avector = [a,bsorted(1:end-1)];
qsub = arrayfun( ...
    @(a,b)integral(f,a,b,'AbsTol',atol,'RelTol',rtol), ...
    avector,bsorted);
qsorted = cumsum(qsub);
q(idx) = qsorted;

But do note that quadgk and integral control the overall error, and doing subintegrations and adding may yield a lower accuracy result than either code would have given in one shot. -- Mike

[Tobias]

arrayfun seems to be what I was looking for. I'll also try the cumsum approach. Thanks a lot.