Trapezoidal integration without built-in functions

Hi guys, I was trying to write a simple code that implements the trapezoidal integration of a f(x,y) function over x and y domain.

I'm wondering if there's anyway to improve the efficiency of that code, to perform it way faster (as the number of points increases exponentially), also because I would like to extend to 3D. In particular I'm wondering if there's any possibility to do it in a single double for and not three for. Thank you so much in advance.

clc; clear all; close all;
x = linspace(-1,1,300);
y = linspace(-1,1,40);
dx = x(2)-x(1);
dy = y(2)-y(1);
[X,Y] = meshgrid(x,y);

f = exp((X-1.5).^2+Y.^2);
fy = zeros(length(y),1);
for j = 1:length(y)
    for i = 1:length(x)-1
fy(j) = fy(j) + (f(j,i)+f(j,i+1))*(x(i+1)-x(i))/2;    
    end
end
int = 0;
for j = 1:length(y)-1
int = int + (fy(j+1)+fy(j))*(y(j+1)-y(j))*0.5;
end

Answer 1

[Torsten]

f = @(x,y) exp((x-1.5).^2+y.^2);
x = linspace(-1,1,300);
y = (linspace(-1,1,40)).';
F = f(x,y);
value_integral1 = trapz(y,trapz(x,F,2))

value_integral1 = 336.7567

value_integral2 = trapz(x,trapz(y,F,1))

value_integral2 = 336.7567

So you see that the one-dimensional trapz can be applied twice. If you organize your code such you will be most effective.

Answer 2

[John D'Errico]

The fact is, trapezoidal integration in a higher number of dimensions in R^d, MUST take O(N^d) operations. Sorry, but it will be an exponential thing, where N is the number of points taken in each dimension, and d is the number of dimensions. This is the nature of the beast.

You can gain some by the use of a tool like ndgrid, so you can avoid needing to write nested loops. And if your code is fully vectorized, you then essentially move the loops to become internal, implicit loops. And that is faster. For example:

v = 0:0.05:1;
[X,Y,Z] = ndgrid(v);
fun = @(x,y,z) x.^2 + y.^2 + z.^2;

% evaluate at all points in one call, so 21^3 = 9261 nodes in total.
W = fun(X,Y,Z);

% Now, you can just call trapz 3 times.
format long g
trapz(v,trapz(v,trapz(v,W,3),2),1)

ans = 
                   1.00125

How well did I do? Analytically, the integral should be 1, so this worked about as well as a trapezoidal rule would be for that set of nodes.

syms Xs Ys Zs
int(int(int(Xs.^2 + Ys.^2 + Zs.^2,Xs,[0,1]),Ys,[0,1]),Zs,[0,1])

ans = 

1

Can you write your own code to perform the integration? Yes, but why? You won't gain much in speed, if anything, and you take the risk of bugs in your own code, maintaining your own code, etc. This is something that has been discussed a hundred times on the forum.

Can you do a better job of the integration? Well, that is likely. For example, you might have chosen to build this in the form of an iterated Gauss-Legendre integral. That would have taken a trivial amount to implement, and it will typically be far more accurate than a trapezoidal integral.

Let me see, I wrote a tool called gaussuadrule long ago (part of my sympoly toolbox.). It tells me the Gauss-Legendre rule for only 3 points has nodes and weights of

Nodes = [ -0.774596669241483  0  0.774596669241483]

Nodes = 1×3

-0.774596669241483 0 0.774596669241483

Weights = [ 0.555555555555556  0.888888888888889  0.555555555555555]

Weights = 1×3

0.555555555555556 0.888888888888889 0.555555555555555

Now, tranform the nodes and weights to the [0,1] unit cube, wher the integration is to be performed.

Nodes01 = (Nodes + 1)/2

Nodes01 = 1×3

0.112701665379259 0.5 0.887298334620741

Weights01 = Weights/2;
[X,Y,Z] = ndgrid(Nodes01);
[Wx,Wy,Wz] = ndgrid(Weights01);

fxyz = fun(X,Y,Z);
I = sum(fxyz.*Wx.*Wy.*Wz,'all')

I = 
         0.999999999999999

It lost just a tiny bit of accuracy, because I pasted in the nodes and weights. So the least significant bits of those values were not truly correct. But even so, the result would be essentially accurate to full precision, since a 3 node Gauss-Legendre rule will be exact for this case. (In fact, it is overkill for this simple problem.)

The point is, this integration required only 27 points in total, instead of 9261.

My point? If you really want to increase the speed, use a better method.