Re: Scattered array, 2D integration?

[Nicholas Jankowski]

On Sat, Oct 5, 2019 at 10:02 PM Nicholas Jankowski <address@hidden> wrote:
>
> > With current functions implemented in Octave, is my only option to create 
> > such a grid a la meshgrid and griddata?
> > Matlab 2013 introduced a scatteredInterpolant function
>

well, potentially solving my own problem, I came across one answer to
this stackexchange:
"How to integrate over a discrete 2D surface in MATLAB?"
https://stackoverflow.com/a/27987120/4598449

In particular, the last part of the answer:

"If your integration area is rather complicated or has holes, you
should consider triangulating your data.

Let's say your integration area is given by the triangulation trep,
which for example could be obtained by trep =
delaunayTriangulation(x(:), y(:)). If you have your values z
corresponding to z(i) = f(trep.Points(i,1), trep.Points(i,2)), you can
use the following integration routine. It computes the exact integral
of the linear interpolant. This is done by evaluating the areas of all
the triangles and then using these areas as weights for the
midpoint(mean)-value on each triangle.

modifying the script a bit (octave doesnt' have delaunayTriangulation
implemented as a class, but it does have delaunay)

given a set of data z(x,y) as three vectors x, y, z:


function int = tri_integ (x,y,z)
  P = [x(:),y(:)]; T = delaunay(x(:), y(:)); z = z(:);
  d21 = P(T(:,2),:)-P(T(:,1),:);
  d31 = P(T(:,3),:)-P(T(:,1),:);
  areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
  int = areas'*mean(z(T),2);
endfunction

testing it out on an irregular array:

r = logspace(-3,0)(2:end-1); theta = linspace(0,2*pi,21)(1:end-1);
x = r.*cos(theta'); y = r.*sin(theta');
[x2,y2]=meshgrid(linspace(-1,1,21);
x =[x(:);x2(:)];y = [y(:);y2(:)];

f = @(x,y) 1-sin(x.^2+y.^2);
z = f(x,y);

tri_integ(x,y,z)
ans =  1.754037103916947
(running 100x = 1.7seconds)

integral2(f,-1,1,-1,1)
ans =  1.754838406712362
(running 100x = 0.6 seconds)

trying the griddata approach:
[x3,y3]=meshgrid(sort(x),sort(y));
z3 = griddata(x,y,z,x3,y3);
trapz(sort(y),trapz(sort(x),z3,2))
ans =  1.753869939077482
(running 100x = 241 seconds)

not much less accurate, but goes from handling >7500 element arrays to
~1.9M element arrays.  the delay comes from the griddata step.  That
would explain hitting memory limits with large meshes.

Anyway, sorry for the message clutter. At least this will be in the
message archive when I need to look it up again. Discussion of
Matlab's scatteredInterpolant mentioned triangulation a few times, so
I'm wondering if it might have done something similar (rather than
just wrapping griddata) for efficiency.  Might revisit that patch once
I'm done this current project.

nickj