Integral2 and array valued functions

I want to calculate the following:

    Og = [4 3;3 4];
    Oginv = inv(Og);

    syms s1 s2
    z = [s1 s2]*Oginv*[s1;s2];

    fun  = matlabFunction(z*exp(-z/2));
    cdf = integral2(fun,-inf,3,-inf,3, 'AbsTol', 1.e-6)

    fun2  = matlabFunction(z^2*exp(-z/2));
    cdf2 = integral2(fun2,-inf,3,-inf,3, 'AbsTol', 1.e-6)

So, in this case I want to evaluate two integrals. For sure I can express it as one integral with (z+z^2)*exp(-z/2), however, my problem is more complicated than this one, so if I do in this way it takes to long to evaluate the integral. I know that I can evaluate integrals for array valued functions with a function integral. Can I do the same with integral2 and apply it to this problem?

I would be very grateful for any help.

Kind regards, Kolibris

Answer 1

[Mike Hosea]

Not in the same way. INTEGRAL2 does not support an 'ArrayValued' option. I prototyped it, but testing revealed that it did not have the performance benefit you might expect. The reason is that you begin to lose the benefits of the adaptive algorithm as each component has a say in where the mesh needs to be refined. The savings from doing the integrals simultaneously are offset by the losses of vectorization within each component and from doing more work in regions where most components do not require it.

If you had a lot of components, I would suggest using ARRAYFUN, or just writing a loop, but with 2 (or several) components you might as well just do as you have done.

Here's another tweak that you might consider. Perhaps you can "tune" it for maximum effect. The 'tiled' method is generally faster than the 'iterated' method in INTEGRAL2, but if you have infinite limits, you have to use the 'iterated' method. So divide your region into parts so that you only need to use the iterated method on the tails.

c = -10;
cdf1 = integral2(fun1,-inf,c,-inf,c)
cdf1 = cdf1 + integral2(fun1,c,3,-inf,c)
cdf1 = cdf1 + integral2(fun1,-inf,c,c,3)
cdf1 = cdf1 + integral2(fun1,c,3,c,3)
cdf2 = integral2(fun2,-inf,c,-inf,c)
cdf2 = cdf2 + integral2(fun2,c,3,-inf,c)
cdf2 = cdf2 + integral2(fun2,-inf,c,c,3)
cdf2 = cdf2 + integral2(fun2,c,3,c,3)

 Naturally you can use a different cut point for y than x if you like.

I'm starting to think we should have errored when people supply only 'AbsTol'. You probably need to supply 'RelTol'. The 'AbsTol' only matters if the result is less than AbsTol/RelTol in magnitude. That is to say, AbsTol is just a safety valve for relative error control when the solution is small in magnitude (since relative error is undefined when the true solution is zero).