2D Simpsons Composite integration
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
In the Mathematics Stack Exchange I questioned about the extension of the composite Simpsons rule to 2D integration.
The problem is that Simpsons rule requires and even number of subintervals and I'm concerned with the possibility of having an even number of subintervals. In another forum, for 1D integration with Simpsons rule, the suggestion was to take out the last
points and use Simpsons 3/8 rule on this last 4 points to keep the integration in the same level of accuracy (the same error order of
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/224175/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/224176/image.png)
My question on the forum was about the extension of this approach to a 2D integration, if the weights were corrected. Assuming they are correct, I coded the function:
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
The question is simple: can I improve this code?
How could I write it without so many if-elseif to account for the different possibilities on the parity of Nx and Ny? How to better define the weights vectors wx and inwy?
0 Commenti
Risposte (0)
Vedere anche
Categorie
Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!