Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Problem with an integral.

Subject: Problem with an integral.

From: Andres

Date: 4 Mar, 2013 21:14:08

Message: 1 of 2

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', 1e-12);
    end


    %DECLARATION OF q2 and Q2


    function y = q2(w, tt)
        if (tt*w) < 1e-5
            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', 1e-12);
    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), 1e-3, t);
    end
    
    Rate(2)

I appreciate any feedback you can give me.

Thanks!

Subject: Problem with an integral.

From: Mike Hosea

Date: 29 Mar, 2013 19:33:08

Message: 2 of 2

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
vector-valued integral with just 1 component. When integrating
vector-valued integrands, INTEGRAL only passes scalars to the integrand
function. To make that work, simply append the option-value 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. 1e-5 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 1e-6, so
when combined with the default RelTol, an AbsTol of 1e-12 only matters
if |Q| < 1e-12/1e-6 = 1e-6.

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 cross-product 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', 1e-12); end
>
> %DECLARATION OF q2 and Q2
>
>
> function y = q2(w, tt)
> if (tt*w) < 1e-5
> 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', 1e-12); 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), 1e-3, t);
> end
> Rate(2)
>
> I appreciate any feedback you can give me.
> Thanks!

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us