MATLAB: Line integral over a scalar field

line integralmatrix indexingscalar field

I have an m by n matrix 'A' full of real values. I need to find the line integral along a line y=mx+b through the matrix. I'm having trouble doing this without for loops. Currently I am creating a zero matrix 'B', then looping through the X indices, solving for the Y index, and setting B(x,y) = 1. I then do an element wise multiplication of A and B and sum the elements of the result. ex. sum(sum(A.*B))
This is a horribly inefficient way of doing things but I can't figure out how to do it quicker. I can easily create arrays of the x,y values that correspond to the line but don't know how to pairwise reference those elements alone.

Best Answer

Um, not only is what you are doing inefficient but also horribly inaccurate, since it forces the solution to lie exactly at some element of that array. So your "line" is actually a jagged mess of an approximation to a line.
Think of the indices of the matrix as floating point values, that interp2 can evaluate. Use interp2, then use trapz. Literally not much more than two efficient lines of code, though setting the problem up will take a couple more lines. And since I'll have MANY comments along the way to explain what I've done...
A = magic(25);
So, x and y are defined over the rectangular domain [1,25]x[1,25] in the (x,y) plane.
I'll arbitrarily pick a line that passes through the square, from one edge to another. Say, y=1.5+2.5*x.
m = 2.5;
b = 1.5;
Where does the line intersect the square? Really, this takes nothing more than testing the 4 edges of the rectangle (here a square) to see where the line intersects.
When x=1, we have y=4. So [1,4] is one end point of the line.
When x=25, we have y=1.5+2.5*25 = 64. This is outside of the square. So we know the line does not pass through the right edge of the domain.
We can see the two possible end points of the line segment as
m*[1,25] + b
ans =
     4    64
So only the first one is relevant here, at (1,4) in the (x,y) plane.
Similarly, when y=1, we solve for x = (y-b)/m. A careful programmer would worry if m could be zero.
([1 25]-b)/m
ans =
       -0.2          9.4
So that locates the point where the line exits the square, as (9.4,25) in the (x,y) plane. When y was 1, the line did not cross the edge, at x=-0.2, so we ignore that case.
I could have located those points programmatically. Regardless, we have identified a line segment with end points at [1,4], and [9.4,25]. It is simple enough to simply generate a relatively fine set of points along the line now, then interpolate using interp2. This won't be exact, because we are only using trapz to do the line integration. You might recognize that along the line we have posed, even a bilinear interpolant does not produce a linear interpolant. It will in fact produce a piecewise quadratic function along that line. (Odds are, you don't care about or know that fact.)
xy1 = [1 4];
xy2 = [9.4 25];
t = linspace(0,1,250)';
xyseg = (1-t)*xy1 + t*xy2;
[xmesh,ymesh] = meshgrid(1:25,1:25);
zseg = interp2(xmesh,ymesh,A,xyseg(:,1),xyseg(:,2));
Now just use trapz on that curve to integrate it. Be careful though. We need to integrate with respect to distance along the curve in the (x,y) plane. Since m could be zero, I'll compute the distance in the (x,y) plane as:
d = cumsum([0;sqrt(sum(diff(xyseg).^2,2))]);
trapz(d,zseg)
ans =
       7608.8
So the line integral as desired is 7608.8. The complete set of computations will have taken a tiny fraction of a second in MATLAB.
Related Question