Back. At least a good drive gives me some time to think. :)
There are at least three if not 4 ways I can think of to solve this numerically. An analytical solution seems something of a stretch though, but let me see where this goes first.
The simplest idea I can think of is to convert to polar coordinates, with the origin at the Weibull mode.
Since this form is separable, we see that the mode happens at the point in (x,y)
(x_m, y_m) = ( ((a-1)/a)^(1/a), ((b-1)/b)^(1/b) )
You can do the work as I just did quickly, or be lazy, and look up the Weibull in Wikipedia, which gives that as the mode.
If we view this in polar form, with that as the origin then we can do a 2-d integration. For any value of the polar angle, use fzero to locate the radial location where z crosses the level of interest (I've called it z0 below.) Then just integrate the function
(z(x_m + r*cos(theta), y_m + r*sin(theta)) - z0) *r*dr*dtheta
where theta runs from 0-2*pi, and r from 0 to the crossing point identified by fzero. No, you cannot use integral2 here easily because the upper limit in r is obtained from fzero, but so what? It is just a nested integration, thus nested calls to integral.
I'd suggest the above algorithm is simplest, as well as reasonably efficient. It won't give you a symbolic solution, but my gut tells me that may not be trivial, if at all possible. And somehow I doubt a symbolic solution for general a & b is even possible, but that is just my gut talking.
My second choice of algorithm seems moderately interesting, a continuation approach. Find one point on the z contour, yielding z0. I'd suggest you again use the polar form I describe above, solving for a single point at theta==0. That won't be that difficult to find. Thus, for y=y_m at the mode in y, use fzero to find where z(x,y) crosses z0.
Now, the entire contour can be derived as a solution to a nonlinear differential equation starting from that point. Thus, given an initial value of r at theta==0, start with this form:
z(x_m + r*cos(theta), y_m + r*sin(theta)) - z0
Differentiate with respect to theta, treating r as a function of theta. Use an ODE solver to walk around the contour line, until you return to theta=2*pi. Now you have the entire contour delineated in the form of r(theta), so the integration becomes almost trivial. Might take a little thought to implement, but it should be doable.
So those are the solution approaches I'd suggest are of interest. Given the contour line as a polygon, in the form of (r(theta),theta), so therefore (x(theta),y*(theta)), then there are many ways to do the integration. You mentioned finite elements, so I assume you understand FEM. Just dissect the region bounded by the convex polygon as a triangulaton, making sure that one vertex of the triangulation lies at the mode point. Now compute the integral using a Gaussian integration with each triangle of the triangulation, compute the integral as the sum of all pieces.
I even have some code that does the Gaussian integration for you of a general function over a triangulated domain, though it is not on the File Exchange. I let it drop off some time ago, but I've attached it (tessquad) to this comment. tessquad does a gaussian integration (of user specified order) within each simplex of a triangulation, then computes the aggregate. It will be fairly accurate, since you can specify any order you want, and this integral will not be too nasty. Tessquad does use another tool from the FEX, simplexquad , by Greg von Winckel. You would need to download it. So really all you need do is build a triangulation from a polygon. Luckily, that too is available. I'd start here with mesh2d, although I did write some meshing and mesh refinement tools of my own. Mesh2d was pretty good. Lots of ways to do this part though.
Best Answer