This question is an extension of this question for 2D integration.
The formulation of the problem is based on this page
Basically, the composite Simpson's rule for 2D integration is
$
\iint_R f(x,y) dx dy\approx\sum_i\sum_jw_{ij}f(x_i,y_j)
$
over the rectangle $R={(x,y):a \le x \le b, c \le y \le d}$.
The domain is subdivided in an even number of intervals with spacing
$h_x=\frac{b-a}{2m},\hphantom{-} x_i=a+ih_x, \hphantom{-} i=0,1,\cdots, 2m \\
h_y=\frac{c-d}{2n}, \hphantom{-} y_j=c+jh_y,\hphantom{-} j=0,1,\cdots,2n$
And the integration is performed using Simpsons rule in each dimension.
For a 1D integration, the composite Simpsons rule can be written as:
$
\int_a^bf(x)dx \approx \frac{h}{3} \left[ f(a)+2\sum_{i=1}^{m-1}f(x_{2i})+4\sum_{i=1}^{m}f(x_{2i-1})+f(b) \right]
$
To implement this, I can create the vector of weights
$
\left[ \matrix{ 1 & 4 & 2 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1 } \right]
$
And calculate (assuming the element-wise multiplication is .*:
I = sum(w.*f)*h/3
For a 2D, the pattern can be extended to a grid in the rectangle:
$
\left[ \matrix{ 1 & 4 & 2 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1 \\
4 & 16 & 8 & 16 & 8 & 16 & 8 & \cdots & 8 & 16 & 4 \\
2 & 8 & 4 & 8 & 4 & 8 & 4 & \cdots & 4 & 8 & 2 \\
\vdots \\
2 & 8 & 4 & 8 & 4 & 8 & 4 & \cdots & 4 & 8 & 2 \\
4 & 16 & 8 & 16 & 8 & 16 & 8 & \cdots & 8 & 16 & 4 \\
1 & 4 & 2 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1 } \right]
$
Therefore, the 2D integrations, given the weights above is:
I = sum(w(:).*f(:))*hx*hy/9
If, however, I have an odd number of intervals, I have to change my functions accordingly. Following this answer, to keep the same order of Simpsons rule (with $m$ odd):
$
\int_a^bf(x)dx, \hphantom{-} h_x=\frac{b-a}{m}, \hphantom{-} x_i=a+ih_x,
\hphantom{-} i=0,1,\cdots,m
$
The idea is to use Simpson's rule for the first $m-3$ points and cover the remaining gridpoints with the Simpsons 3/8 rule. With the weight vectors defined above, we have
w_simpson=ones(1,m-3);
w_simpson(2:2:end-1) = 4;
w_simpson(3:2:end-2) = 2;
w_38 = [1 3 3 1];
I = sum(w_simpson.*f(1:m-3))*h/3 + sum(w_38.*f(m-3:m))*3*h/8
Now, the question: to expand this 3/8 rule for a 2D integration. Are the matrices defined below correct?
$
\left[ \matrix{ 1 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1 \\
4 & 16 & 8 & 16 & 8 & \cdots & 8 & 16 & 4 \\
2 & 8 & 4 & 8 & 4 & \cdots & 4 & 8 & 2 \\
\vdots \\
2 & 8 & 4 & 8 & 4 & \cdots & 4 & 8 & 2 \\
4 & 16 & 8 & 16 & 8 & \cdots & 8 & 16 & 4 \\
1 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1 } \right]
\left[ \matrix{3&3&1\\
12&12&4\\
6&6&2\\
\vdots\\
6&6&2\\
12&12&4\\
3&3&1
} \right]\\
\left[ \matrix{3 & 12 & 6 & 12 & 6 & \cdots & 6 & 12 & 3\\
3 & 12 & 6 & 12 & 6 & \cdots & 6 & 12 & 3\\
1 & 4 & 2 & 4 & 2 & \cdots & 2 & 4 & 1\\} \right]
\left[ \matrix{9&9&3\\9&9&3\\3&3&1} \right]
$
If they are correct, then the 2D integral is:
x = a:hx:b;
y = c:hy:d;
[X,Y] = meshgrid(x,y);
Z = fun(X,Y);
I = sum(sum(w1.*Z(1:n-3,1:m-3)))*hx*hy/9 + sum(sum(w2.*Z(1:n-3,m-3:m)))*hx*hy/8 + ...
sum(sum(w3.*Z(n-3:n,1:m-3)))*hx*hy/8 + sum(sum(w4.*Z(n-3:n,m-3:m)))*hx*hy*9/64
Is it correct? If so, is this the best way to implement the integration?
How do I define the weight matrices?
In Matlab, we must use the sum function twice, because I set a range on matrix. Is there a better way to code it without using for loops?
Edit:
The code came out as this:
function I = simpson2D(x,y,z)
% Composite Simpson's rule for 2D integration
from meshgrid
Nx = size(x,2); % number of columns
Ny = size(x,1); % number of rows
hx = x(1,2)-x(1,1);
hy = y(2,1)-y(1,1);
if (mod(Nx,2)) % Nx is odd (even number of intervals)
NxIsEven = false;
nx = Nx;
wx = ones(1,Nx);
wx(2:2:nx-1) = 4;
wx(3:2:nx-2) = 2;
else % Nx is even (odd number of intervals)
NxIsEven = true;
nx = Nx-3;
wx = ones(1,Nx);
wx(2:2:nx-1) = 4;
wx(3:2:nx-2) = 2;
wx(Nx-2:Nx) = [3 3 1];
end
if (mod(Ny,2)) % Ny is odd (even number of intervals)
NyIsEven = false;
ny = Ny;
wy = ones(Ny,1);
wy(2:2:ny-1) = 4;
wy(3:2:ny-2) = 2;
else % Ny is even (odd number of intervals)
NyIsEven = true;
ny = Ny-3;
wy = ones(Ny,1);
wy(2:2:ny-1) = 4;
wy(3:2:ny-2) = 2;
wy(Ny-2:Ny) = [3 3 1];
end
w = wy*wx;
if (NxIsEven && NyIsEven)
I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
(sum(sum(w(1:ny,nx:Nx).*z(1:ny,nx:Nx)))+sum(sum(w(ny:Ny,1:nx).*z(ny:Ny,1:nx))))/8 + ...
sum(sum(w(ny:Ny,nx:Nx).*z(ny:Ny,nx:Nx)))*9/64;
elseif (~NxIsEven && NyIsEven)
I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
sum(sum(w(ny:Ny,1:nx).*z(ny:Ny,1:nx)))/8;
elseif (NxIsEven && ~NyIsEven)
I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
sum(sum(w(1:ny,nx:Nx).*z(1:ny,nx:Nx)))/8;
else % (~NxIsEven && ~NyIsEven)
I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9;
end
I = I*hx*hy;
end
Any hints on how to improve it (to treat the parity of Nx/Ny without so many if-elseif-else and code repeating on the definition of the weight matrices wx and wy?
Best Answer
It can be written with some complicated expressions, but I don't know whether you think it would be an improvement:
The idea is to make the Simpson $1/3$ pattern with an extra zero appended if there is an even number of data points. Then we add an equal-length array which is all zeros for an odd number of data points but corrects the last three Simpson $1/3$ rule and the appended zero to the Simpson $3/8$ rule for the last four data points for an even number of data points. Investing a lot of computational effort on getting
wxandwyshould not be a problem because they are $O(N)$ processes while the $2$-d integral is $O(N^2)$.I think your matrices are incorrect because you left out the first point of the Simpson $3/8$ rule parts, but your code right under your matrices looks OK at first glance, although I think Matlab should reject your final block of code. Does it?
EDIT: I found a way to make my answer perhaps a little more palatable. I created two helper functions, the first of which creates the Simpson'2 $1/3$ weights using Matlab's answer to an ac-implied-do and the second uses a Matlab replacement for the ternary operator to choose between an extension of the weights by Simpson's $3/8$ rule. Then the functions are invoked with the proper arguments to produce
wxandwy.The code, especially for helper function
fis still pretty dense but did meet the requirement of eliminatingforloops and replacedif-elseif-elsesyntax with arithmetic expressions likefloor((N-3)/2)andmod(N,2). Notice that the $2$-d Simpson's rule can be expressed as a contraction of the array of functional values so there is no need to form a $2$-d array of weights.