Scattered array, 2D integration?

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Scattered array, 2D integration?

nrjank
I have model data that is a non uniform mesh over a rectangular area (x,y), and a single nodal value at each point. (The data is a n x 3 point-value array, [x,y,S]. I need to do a 2D numerical integration over that domain.  If it was a full rectangular grid, even if nonuniform, I could just call trapz twice.  With current functions implemented in Octave, is my only option to create such a grid a la meshgrid and griddata? Level of refinement on one corner of the mesh makes interpolating to a uniform mesh as per the scattered data example in the manual [1] fairly unwieldy (mapped mesh without loss of refinement can square the amount of data in each vector).

Matlab 2013 introduced a scatteredInterpolant function that they recommend using as an intermediary between the data and other integration functions, but no compatible option implemented in octave yet.  

Was hoping someone would have a less resource intensive method, or a good manual implementation of the interpolant scheme.


Nick J.


Reply | Threaded
Open this post in threaded view
|

Re: Scattered array, 2D integration?

nrjank
> 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

just came across bug #35821 (https://savannah.gnu.org/bugs/?35821)
which has a ?unreviewed? patch from 2013 that is a wrapper for
griddata.  So far my manual tests with griddata are running into
memory issues since it forces a full rectangular mesh from my sparse
data.  I assume I'd run into the same issue with this wrapper.  Is
anyone familiar with the scatteredInterpolant (or TriScatteredInterp)
functions in Matlab?  Does memory use make it seem like they are using
the sparse data or hiding a full rectangular array under the hood?


Reply | Threaded
Open this post in threaded view
|

Re: Scattered array, 2D integration?

nrjank
On Sat, Oct 5, 2019 at 10:02 PM Nicholas Jankowski <[hidden email]> 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