2D quadrature for an array valued function

I need to perform integration of an array valued function over a rectangular domain. I have found dblquad and quad2, which will perform 2d integration with a scalar function, and quadv, which will perform 1d integration of an array valued function, but I need to do both. I put together a couple for loops that will perform a crude rectangular integration but it is quite slow. I was hoping to optimize my code with the built-in functions. Is there any way to combine the above functions to achieve what I'm looking for?

.

edit: Another approach would be to redefine the array valued function as a set of scalar functions, integrate them separately using dblquad, then reassemble them into an array. I'm working with a 12 x 12 array, however, which would mean defining 144 separate functions. Is there an efficient way to do this without writing out each function?

[Doug Hull]

What does it mean to integrate over an array valued function over a rectangular domain? So at each x,y coordinate you would have an (n x m) result?

[Sam G.]

Yes, if you evaluate the function at a position (x,y), you would get an m x n array in return. In addition, if you integrate an array valued function over a rectangular domain, you get an array in return.

[Doug Hull]

What does it mean to integrate in this context? If this returned a scalar at each (x,y), I would integrate the scalar. How do I integrate an m x n array? Would you treat this as m*n different scalars to integrate separatly?

Best practice is to edit your question to reflect this new information. Then I will delete my comments for clarification, and you can too.

Answer 1

[Sam G.]

Well, if anyone is interested in this...

After some more looking around, I found an example of some code that I could adapt to my application. It's basically just 2D gaussian quadrature but since it's my own code I can make it compatible with array valued functions.

% Define weights and abscissas for guassian quadrature
n = 3;
switch n
    case 2
        z = [-0.57735026918962576451 0.57735026918962576451];
        w = [1 1];
    case 3
        z = [-0.77459666924148337704 0 0.77459666924148337704];
        w = [0.55555555555555555556 0.8888888888888889 0.55555555555555555556];
    case 4
        z = [-0.86113631159405257522 -0.33998104358485626480 ...
            0.33998104358485626480 0.86113631159405257522];
        w = [0.34785484513745385737 0.65214515486254614263 ...
            0.65214515486254614263 0.34785484513745385737];
    case 5
        z = [-0.90617984593866399280 -0.53846931010568309104 0 ...
            0.53846931010568309104 0.90617984593866399280];
        w = [0.23692688505618908751 0.47862867049936646804 0.56888888888889 ...
            0.47862867049936646804 0.23692688505618908751];
end

% Defined domain of integration in y
a = 10;
b = 20;

% Define domain of integration in x
c = 50;
d = 55;

% Initialize solution matrix
solution = zeros(3,3);

% Integration in y direction
for i = 1:n
    % Translate abscissa to from [-1 1] domain to [a b] domain
    y = (a + b)/2 + (b -a)/2*z(i);
    % Initialize matrix for the x interval
    xint = zeros(3,3);
    % Integration in x direction
    for j = 1:n
        % Translate abscissa from [-1 1] domain to [c d] domain
        x = (c + d)/2 + (d - c)/2*z(j);
        % Evaluate the function
        fxy = f(x,y);
        % Summation of the weighted function values
        xint = xint + w(j).*fxy;
    end
    % Translate to [-1 1] domain for summation
    xint = (d - c)/2.*xint;
    % Summation of the weighted function values
    solution = solution+ w(i)*xint;
end
% Translate to [-1 1] domain for summation
solution= (b - a)/2.*solution;

Where

% Define the function to integrate
function Q = f(x,y)

Q(1,:) = [-6.*x -2.*y -6.*x.*y];
Q(2,:) = [-2.*x -6.*y -6.*x.*y];
Q(3,:) = [-4.*y -6.*x^2 -6.*y^2];

[Sam G.]

As you can see, the example function returns an array and the result of the integration is also an array. This type of calculation is useful in 2D finite element solvers, which is my application (this code has been generalized a bit for easier understanding). I still don't get the advantage of pre-compiled code but it is an order of magnitude faster without any loss of accuracy.