MATLAB: How to determine the volume bounded by z=x^(a-1)y^(b-1)exp(-x^a-y^b) and z=0.05 (a,b>1) x>=0 y>=0

contourintegrationnumerical integration

I have difficulty computing the volume of the above problem. As the function is implicit, I am unable to write the boundaries in terms of function.

I have used integral2, integral3, dblquad and other functions but they obviously do not give correct answers as limits must be functions.

The volume may be numerically computed. Is there a script that lets you enter these equations and computes the volume itself? If not, how do I calculate it?

Best Answer

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.

Related Solutions

MATLAB: How to plot a 3D surface bounded over the first octant

It fails to do what you asked to do, because what you did has nothing to do with constraining it to the first octant. You generated a list of (x,y) pairs, but provided no understanding to MATLAB of the real question, which is what is the surface in terms of where thing lives in the first octant!

If you recognize the result will be a planar triangle, embedded in R^3, then you could simply compute the coordinates of the three vertices of that triangle. Then use patch.

So when x=y=0, what is z? z=a of course. Likewise, when x=z=0, what is y? Finally, when y=z=0, what is x?

You should see that this results in three vertices of a triangle, thus the rows of

[0 0 a;0 a 0;a 0 0]

Patch is now simple to use.

help patch

If your goal is somewhat less understanding of the result, so you hope that MATLAB can simply be told to plot some general relation when constrained to the first octant, that will take a completely different approach, because then you are effectively asking MATLAB to do your thinking for you.

MATLAB: The length of the intersection curve between a sphere and a cylinder

I won't do your homework. I'll get you off to a good start, in fact, probably more than you deserve, since you appear to have made no serious effort.

The sphere is a sphere centered at the origin, of radius 4.

The cylinder is a circle in the (x,y) plane, so the same for all z. Completing the square, it can be written as:

(x^2 - 2*x +1) + y^2 = 1

(x - 1)^2 + y^2 = 1

So the cylinder has center (1,0,z), for all z.

That tells me the two objects do intersect, as well, they will intersect in TWO disjoint curves, one for positive z, the second for negative z.

In fact, if you view the intersection curve from the top, it will look perfectly circular. And this should be the clue you need.

Suppose you thought about this in terms of polar coordinates? To make it simple, translate the problem so the cylinder is centered at (0,0) in (x,y). Then the equation of the sphere is just:

(x+1)^2 + y^2 + z^2 == 16

The cylinder is now

x^2 + y^2 == 1

Nothing changes about the length of the curve. But perhaps now you can visualize that happens when we re-write this all in polar form. Actually, cylindrical coordiantes, thus (r,theta,z)

In cylindrical coordinates, the cylinder has equation

r == 1

Now, substitute into the sphere equation. Remember, since we now know that r==1,

x = r*cos(theta) = cos(theta)

y = r*sin(theta) = sin(theta)

now, lets look at what happened to the sphere.

(I've deleted the stuff that gives you a formula for z as a function of theta. In fact there is a remarkably simple form for z that will result.)

z = stufffffffffffff

Can you now find the arc length of that function, as theta varies from 0 to 2*pi? (I hope so.) A quick numerical computation of my own had the arc length as roughly 6.3864. My guess is only the first few digits of that are accurate, because of the approach I took for that first pass. But it should be close to 2*pi, but a bit larger, just as a visual assessment from a plot.

How do you find the arc length of a curve? I might be lazy, and use my arclength utility, as posted on the file exchange.

theta = linspace(0,2*pi,10000);
x = cos(theta);y = sin(theta);z = ????????;
arclength(x,y,z,'spline')
ans =
          6.39448891369661

You should be able to do the same however, using an integration. A google search for the phrase

"arc length of a parametric curve"

might give you some ideas. Don't forget that this curve is in the form (x(theta),y(theta),z(theta)). x(theta) is just cos(theta), etc. In fact, when I do that computation, integral tells me the arc length should be

ans =
          6.39448891369662

int did not find a symbolic solution however. Que sera, sera. But vpaintegral gives me:

ans =
6.3944889136966193144

I would note that the answer is indeed just a wee bit more than 2*pi, as I predicted it should be.

So, now all you need to do is do the algebra, then do some numerical methods work to replicate my results. IF you do decide to make a serious effort, then show what you did as a comment here. If you come at least close, showing that you have thought about the problem, then I might raise some of the veil on my answer.

Best Answer

Related Solutions

MATLAB: How to plot a 3D surface bounded over the first octant

MATLAB: The length of the intersection curve between a sphere and a cylinder

Related Question