Numerical integration using Trapz

Hi all,

I am trying to calculate a double integral over a rectangular grid using trapz.

The grid is defined as:

[x,z] = meshgrid(xp,zp);

xp and zp are the 1d arrays of x and z coordinates of the grid.

I trying to calculate the function values (F) at the grid nodes as:

F(i,j)=trapz(zp,trapz(xp,g.*f,2));

here,

g is a function of x and z i.e. g(x,z) whereas, f values are already calculated values at each node of the grid i.e. not a function of x and z.

Can someone suggest how may I do it, please.

Answer 1

[Star Strider]

That should produce one scalar value as ‘F’, so the subscripts are not necessary.

I would do something like this —

F = trapz(zp,trapz(xp,g.*f));

The result of the first trapz call will be a vector, and of the second will be a scalar. So long as ‘g’ and ‘f’ have compatible dimensions, and they match ‘xp’ and ‘zp’ appropriately, that should work.

[Waseem Akhtar]

Thank you for the quick reply.

Actually, I am trying to figure out how may the trapz can take different values of 'f' corresponding to each node/grid point. Unlike 'g', these are not functions of x and z but are already calculated unique values at each grid point.

[Star Strider]

My pleasure!

I am not certain what you want to do. If you want to experiment with different combinations of the elements of ’f’, one way would be to create a separate matrix (perhaps ‘fm’) of ones, then substitute specific elements of ‘f’ into it, then do the integrations.

So for example —

F = zeros(size(f));                                             % Preallocate
for i = 1:size(f,1)
    for j = 1:size(f,2)
        fm = ones(size(f));
        fm(i,j) = f(i,j);
        F(i,j)=trapz(zp,trapz(xp,g.*fm));
    end
end

Assuming that I understand that this is what you want to do, that is likely the most efficient way I can think of to do it.

[Waseem Akhtar]

Thanks once again!

I think we are close, let me try to explain it a bit more.

Consider that I have a grid of size 10x5 (say xp=1:1:10 and zp=1:1:5). And, [x,z]=meshgrid(xp,zp);

I want to calculate the value of 'F' say F(i,j) at each of the 50 nodes such that, for each F(i,j) value, the expression F=trapz(zp,trapz(xp,g.*f)); should return a scalar.

here;

'g' is a function of (x,z) defined as say, g = x.^2 + z.^2;

'f' is the matrix containing already calculated unique values corresponding to each node say, ones(size(f)).

what I want is that each value of 'g' (at each node corresponding to each x and z pair) should get multiplied by the corresponding value of 'f' at the corresponding node during evaluation of F=trapz(zp,trapz(xp,g.*f));

Note that I am using another outer nested loop to get the F(i,j) value at each of the 50 nodes such that i=10 and j=5.

Please bear with my poor explanation, I am not an expert of matlab.

Thank you for your time.

[Star Strider]

As I understand that, my code does exactly what you want. I do not understand the reason for the loops, however, unless it is to do exactly what my code does.

My code multiplies the ‘fm’ matrix by ‘g’, with corresponding values of ‘fm’ first being replaced by selected values of ‘f’. The multiplication is element-wise, so only the element of ‘g’ that is multiplied by the specific value of ‘f’ selected in the loop changes. That result is then reflected in the corresponding element of ‘F’.

This is the best I can do. If it is still not ‘correct’, please provide an example illustrating the necessary calculation to get the desired result.

[Waseem Akhtar]

Thank you Star Strider, your first answer was correct.

[Star Strider]

As always, my pleasure!