For efficiency (and it *really* matters for that), MATLAB tries to
vectorize the calls to the integrand function. If you read the doc, it
will tell you that your integrand function must accept an *array* of
inputs and return an array of corresponding outputs. If you have a
function which cannot do that, and most of your functions cannot, you
have a couple of options. The simplest is to regard your problem as a
vectorvalued integral with just 1 component. When integrating
vectorvalued integrands, INTEGRAL only passes scalars to the integrand
function. To make that work, simply append the optionvalue pair
'ArrayValued',true
to the argument list for each call to INTEGRAL.
The other approach (and almost always a faster one) is to use ARRAYFUN
to vectorize a function that only accepts scalars. So, for example, to
vectorize q2, you might use
q2v = @(w,tt)arrayfun(@q2,w,tt);
instead of q2 in the definition of Q2. I do think ARRAYFUN can get
confusing, so if I were you I would use the 'ArrayValued' option *first*
and then try to remove those one at a time by introducing uses of ARRAYFUN.
I notice that you are only supplying 'AbsTol' below. That is probably a
mistake. INTEGRAL attempts to meet a *mixed* error control condition
involving an estimate of the relative error. Absolute error control is
only relevant when the true answer approaches zero. A good way of
thinking about it is that the 'RelTol' indicates how many significant
digits you want (e.g. 1e5 means about 5 significant digits, or almost
that), and 'AbsTol' is a kind of "threshold", a "good enough" error when
it becomes difficult to achieve the requested relative tolerance because
cancellation of positive and negative areas lead to a small correct
answer. Technically, only 'AbsTol' only matters if Q < AbsTol/RelTol
and only 'RelTol' matters otherwise. The default RelTol is 1e6, so
when combined with the default RelTol, an AbsTol of 1e12 only matters
if Q < 1e12/1e6 = 1e6.
Finally, you will probably observe that when nesting integrals you will
get more accuracy than you ask for. That sounds like a good thing, but
it is not coming for free. There are a couple of reasons for it, one
being that INTEGRAL has a conservative initial mesh (meaning that there
are quite a few evaluations points at a minimum). The other reason I
suspect is more technical and has to do with truncation error associated
with crossproduct terms.

Mike
Andres wrote:
> Hi,
> I am new to MATLAB and I am encountering a bit of trouble writing a
> program to calculate a an integral which itself depends on functions
> that get integrated each time the integrand is called.
> I have already written a function for the integrand, which, when called
> for a given value of its argument, yields the result that I expect.
> However, when I try to integrate over this argument, I get multiple
> errors about disagreeing matrix dimensions. The strange thing is that I
> can evaluate the integrand at all values for the argument.
> My code is as follows:
> %DECLARATION OF q1 and Q1
>
> function y = q1(w, tt)
> y = exp(tau*w).*sin(tt.*w)./w;
> end
>
> function y = Q1(tt)
> y = integral(@(w) q1(w, tt), 0, Inf, 'AbsTol', 1e12); end
>
> %DECLARATION OF q2 and Q2
>
>
> function y = q2(w, tt)
> if (tt*w) < 1e5
> y = (1/24.0)*exp(tau.*w).*tt.*tt.*(12.0 
> w.*w.*tt.*tt).*coth(0.5*w).*w;
> else
> y = exp(tau.*w).*(1.0  cos(tt.*w)).*coth(0.5*w)./w;
> end
> end
>
> function y = Q2(tt)
> y = integral(@(w) q2(w, tt), 0, Inf, 'AbsTol', 1e12); end
>
> %RATE INTEGRAL
> function y = Integrand(tt)
> y = cos(Omega*tt + Q1(tt)).*exp(Q2(tt));
> end
>
> function y = Rate(t)
> y = integral(@(tt) Integrand(tt), 1e3, t);
> end
> Rate(2)
>
> I appreciate any feedback you can give me.
> Thanks!
