2D quadrature for an array valued function

28 views (last 30 days)
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
Doug Hull on 24 Jan 2011
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.

Sign in to comment.

Accepted Answer

Sam G.
Sam G. on 21 Feb 2011
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];
% 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;
% Translate to [-1 1] domain for summation
xint = (d - c)/2.*xint;
% Summation of the weighted function values
solution = solution+ w(i)*xint;
% Translate to [-1 1] domain for summation
solution= (b - a)/2.*solution;
% 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];
  1 Comment
Sam G.
Sam G. on 21 Feb 2011
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.

Sign in to comment.

More Answers (0)


Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!