Composite 2D Simpsons Rule with odd intervals

Question

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?

Answer 1

It can be written with some complicated expressions, but I don't know whether you think it would be an improvement:

wx = hx/3*[1 4 reshape([2;4]*ones(1,floor((Nx-3)/2)), ...
    1,[]) 1 zeros(1,mod(Nx+1,2))]+...
    mod(Nx+1,2)*hx/24*[zeros(1,2*floor((Nx-3)/2)) ...
    1 -5 19 9*ones(1,mod(Nx+1,2))];
wy = hy/3*[1 4 reshape([2;4]*ones(1,floor((Ny-3)/2)), ...
    1,[]) 1 zeros(1,mod(Ny+1,2))]+...
    mod(Ny+1,2)*hy/24*[zeros(1,2*floor((Ny-3)/2)) ...
    1 -5 19 9*ones(1,mod(Ny+1,2))];
I = wy*z*wx'

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 wx and wy should 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 wx and wy.

g = @(N,h) h/3*[1 4 reshape([2;4]*ones(1,floor((N-3)/2)),1,[]) 1];
f = @(N,w) cell2mat(getfield({[w(1:end-3) [w(end-2:end) 0]+ ...
    w(1)/8*[1 -5 19 9]],w},{1+mod(N,2)}));
wx = f(Nx,g(Nx,hx));
wy = f(Ny,g(Ny,hy));
I = wy*z*wx';

The code, especially for helper function f is still pretty dense but did meet the requirement of eliminating for loops and replaced if-elseif-else syntax with arithmetic expressions like floor((N-3)/2) and mod(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.