2D quadrature for an array valued function

11 次查看（过去 30 天）
Sam G.
Sam G. 2011-1-21
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?
  3 个评论
显示 1更早的评论
Sam G.
Sam G. 2011-1-24
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
Doug Hull 2011-1-24
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sam G.
Sam G. 2011-2-21
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];
  1 个评论
Sam G.
Sam G. 2011-2-21
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by