Calculating the height of a liquid in a sphere using secant method

The problem is to determine how much the spherical tank must be filled to contain 30m3 if

R = 3m. The volume is calculated according to equation: V = (𝜋ℎ^2*[3𝑅 − ℎ])/3

Use the secant method. Print the approximations, the relative error. Ensure that the calculation ends when the error is less than a certain value.

Find the correct solutions using the appropriate matlab function. Print approximations and errors.

Bonus: show graphically the dependance of the Volume V on the Height h. f(x).

[James Tursa]

What have you done so far? What specific problems are you having with your code?

[Justin Michaellmore]

This is what I did so far:

clear all
clc
er=100;
es=0.001;
func = @(h)(pi*h^2*(9-h))/3;
h1=1;
h2=2;
f1=func(h1);
counter=0;
while abs(er) > es
  counter=counter+1;
  f2=func(h2);
  er=((h2-h1)*f2)/(f2-f1);
  fprintf('error:%f\n',er); 
  h1=h2;
  f1=f2;
  h2=h2-er;
  fprintf('height (aprox):%f\n',h2);
endwhile

I don't know how to calculate height. When I put h2 into equation, it gets nowhere near 30

Also I don't know what matlab function to use to calculate exact number.

Would appreciate help.

Thank you.

Answer 1

[James Tursa]

You need to use a function that is 0 at the h you are looking for. So you need that 30 in your function handle. E.g.,

func = @(h)(pi*h^2*(9-h))/3 - 30

If you want to compare to a MATLAB function answer, realize that you have a polynomial in h. So just make a vector of the coefficients and feed it to the roots( ) function.

[Justin Michaellmore]

I edited my mistake (kept dividing by 30 instead of subtracting (it's late)), but now I get endless loop. What did I do wrong? Also am I printing correct variables?

[James Tursa]

When I run your exact code above with the -30 appended to your function handle, here is the output I get:

error:-0.032394

height (aprox):2.032394

error:0.005507

height (aprox):2.026888

error:-0.000018

height (aprox):2.026906

And when I plug in the final h2 value:

>> func(h2)

ans =

-3.0165e-07

So things seem to be working for me OK.

[Justin Michaellmore]

Didn't work for me so I ended up using slightly different code:

clear all

clc

er=100;

es=0.001;

func = @(h)((pi*h^2*(9-h))/3)-30;

x1=1;

x2=2;

f1=func(x2);

max=100;

n=2;

h1(n)=x1;

h1(n-1)=x2;

while abs(er) > es

f2=func(x1);

h1(n+1)=((x2*f2)-(x1*f1))/(f2-f1);

ea=abs(h1(n+1)-h1(n));

er=100*ea/h1(n);

fprintf('error:%f\n',er);

fprintf('height (aprox):%f\n',h1);

r=fzero(func,x1);

fprintf('matlab solution: %f\n',r);

rr=h1(n)-r;

fprintf('discrepancy: %f\n',rr);

if n>max

break

endif

n=n+1;

endwhile

Still don't know how to show correlation between volume and height graphically though. This is what I did:

ezplot(func, [0,x1]);

hold on;

y=func(h1);

plot(h1,y,);

Thanks for your help nevertheless

[James Tursa]

"Didn't work for me"

So, that doesn't tell me much. Running your exact code with the -30 addition produced the answer for me. What is it that didn't work for you? You should be getting the exact same answer I did.

clear all
clc
er=100;
es=0.001;
func = @(h)(pi*h^2*(9-h))/3 - 30;
h1=1;
h2=2;
f1=func(h1);
counter=0;
while abs(er) > es
  counter=counter+1;
  f2=func(h2);
  er=((h2-h1)*f2)/(f2-f1);
  fprintf('error:%f\n',er); 
  h1=h2;
  f1=f2;
  h2=h2-er;
  fprintf('height (aprox):%f\n',h2);
end

As for the plot, just realize that the height can range from 0 to 2*R. So create that range and plug it into your vectorized V formula:

R = 3;
h = 0:0.01:2*R;
V = (pi*h.^2.*(9-h))/3;

Then plot(h,V)