Volume of a Curve Revolved about y-axis

Hi all,

I need to find a value b that satisfies

Where r is a known constant and b is unknown. Instead of evaluating this mathematically, I want to use a numerical approximation in MatLab

This is supposed to be the volume created by rotating a quarter of an ellipse about the y-axis

the original equation is

Here is my current code:

clc
clear
%% Calculation of R

mA = 5.8;
global r;
r = sqrt(mA/pi);


%% Numerical Calculation of b
% This section uses the given map area and average elevation to approximate
% a value for b

expectedVol = .0024*5.8;

x = 0:.001:r;
bv = 250;

calcVol = 1;

while calcVol >= expectedVol
    
    y = x.*sqrt(r^2-x.^2)./bv;
    
    area = trapz(x,y);
    
    calcVol = area*2*pi;
    
    bv = bv+.01;
    disp(calcVol);
end

%% Verification of b
% If the above approximation of be is be believed, then the integral to
% evaluate the avergae elevation of the island should be equal to .0024 km

x = linspace(-r,r,1500);
y = sqrt(r^2-x.^2)./bv;

area = trapz(x,y);

havg = area/(2*r);

figure(1)
hold on
    axis([-r*1.1,r*1.1,0,.01])
    pbaspect([2,1,1])
    plot(x,y)
    scatter(0,max(y),'rx')
    plot([-2*r,2*r],[havg,havg],'k')
    legend('Island','Max Elevation','Average Elevation')
hold off

This code, in theory, should start with a value of b = 0, then increase it until the Calculated Volume equals the Expected Volume. The third section checks the value of b by calculating the average value of the function, which should yield .0024, but it doesn't. Am I calculating theese integrals correctly? I know there's an integral function, but I don't know how to apply that with volume.

Thanks for your help,

Ben

[Torsten]

You can directly solve for b - why do you iterate ?

b = 2*pi/0.0139*integral_{x=0}^{x=r}x*sqrt(r^2-x^2) dx

[Benjamin Clements]

Wow, I guess I was so caught up in the problem, I forgot how to do basic calculus. Now the other question is why isn't the average elevation calculation not returning the expected value of .0024?

Answer 1

[John D'Errico]

Solving for b numerically as you want to do is a really bad way to write this. If b enters in in a way that is more complex, then use fzero.

Sorry, but it is. And using a global variable there is as silly, for no good reason. Don't use globals. Ceraintly don't use globals when you have no need for them. And you never have a need for them.

Next, it is virtually never a good idea to solve a problem by stepping through all possible values until you find something that comes close to a solution. Again, at worst, you might use fzero. However...

b is an unknown constant.

syms x r b
int(x*sqrt(r^2 - x^2),x,[0,r])
ans =
r^3/3

So we have

2*pi*r^3/(3*b) = 0.0139

and basic algebra yields

b = 2*pi*r^3/(3*0.0139)

There are many other issues/problems with your code. For example, it is a BAD idea to use construct like this:

x = 0:.001:r;

Why? Because r is not an integer multiple of 0.001.

r = sqrt(mA/pi)
r =
       1.3587

So the last step of that vector will stop at 1.358, and your integral will be less accurate than necessary. Instead, use a construct like this:

nx = 1 + ceil(r/0.001);
x = linspace(0,r,nx);

That will give you a spacing that is approximately 0.001, and that will go all the way up to r.

[John D'Errico]

Did nobody yet answer that followup question I just saw in the comments?

You want to now compute the average "elevation", of a circularly symmetric volume. The volume is given as 0.0139. So the value of b here is irrelevant. All we need to know is the radius, r.

r
r =
       1.3587

0.0139/(pi*r^2)
ans =
    0.0023966

That is, if you have such a volume, the volume is related to the average height as

volume = 0.0139 = area*averageHeight

where area is the area in the plane of that volume, and averageHeight is your goal to compute. (This idea is valid for any such volume.)

The area of a circle is just pi*r^2. So it is 0.0023966..., which seems to be the result you were expecting.

[Benjamin Clements]

Yes, that is the value I'm looking for, but I suppose I failed to mention that that's how the value of .0024 was obtained, I was hoping validate the value of b by using a different method to calculate the average 'elevation'.

[Srikant Upadhyay]

hey bro i have a question too i have posted it https://in.mathworks.com/matlabcentral/answers/1575853-i-am-trying-to-get-volume-of-curve-rotated-but-getting-error-using-fill-and-error-in-volume-curve-li?s_tid=srchtitle could you check it