Bisection method code - question.

Hi, I wrote the following function for solving V=L[arccos(h/r)r^2 - h(r^2-h^2)^0.5] using the bisection method. However, as I execute the program it gets stuck, yet I cannot figure out why. I'd appreciate any comments.

function h = volume(V,L,r,h1,h2)
h=(h1+h2)/2;
err=abs(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2)));
while (err > 0.01)
    if (V-L*(acosd(h1/r)*r^2-h1*sqrt(r^2-h1^2))*(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2))) < 0)
        h2=h;
    else
        h1=h;
    end
    h=(h1+h2)/2;
    err=abs(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2)));
end

Answer 1

[Matt J]

I think you're just missing some parentheses in

 if (V-L*(acosd(h1/r)*r^2-h1*sqrt(r^2-h1^2))*(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2))) < 0)

[Matt J]

It would be easier to read and debug if you calculate "err" using a function instead of repeating its expression all the time,

 function h = volume(V,L,r,h1,h2)

   err=@(h) abs(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2))); %anonymous func

      h=(h1+h2)/2;
      while (err(h) > 0.01)
          if err(h1)*err(h) < 0
              h2=h;
          else
              h1=h;
          end
          h=(h1+h2)/2;
      end

 end

[Matt J]

It's also possible that there is no root err(h)=0 in the interval [h1,h2].

Your code does nothing to check, for example that err(h1)*err(h2)<0 for the initial h1, h2.

[Yuval]

Okay, now the following code doesn't get stuck but the answer seems to be inaccurate. Any idea why?

function h = volume(V,L,r,h1,h2)
h=(h1+h2)/2;
err=abs(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2)));
while (err > 0.01)
    if ((V-L*(acosd(h1/r)*r^2-h1*sqrt(r^2-h1^2)))*(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2))) < 0)
        h2=h;
    else
        h1=h;
    end
    h=(h1+h2)/2;
    err=abs(V-L*(acosd(h/r)*r^2-h*sqrt(r^2-h^2)));
end

I mean, it keeps giving out h=0.9998, even when I change the while condition to (err > 2.009032E-11), for instance.

[Matt J]

What data are you using for V,L,r,h1,h2?