Error, cannot determine if this expression is true or false (Bisection method)

Viewed 289

I am writing code for bisection method, and I keep getting this error with almost every method I use. I can't seem to identify the problem, I am supposed to use bisection method, false method position, and secant method. Everytime I run the code it gives me this error code

                                   -8
                          E := 1 10  

S := 0.1*10^(-7);
                                   -8
                          S := 1 10  

f := x -> 950*ln(200000/(200000 - 3000*x)) + (-1)*9.81*x - 500;
f := proc (x) options operator, arrow; 950*ln(200000/(200000-300\

  0*x))-9.81*x-500 end proc


a := 10;
                            a := 10

b := 50;
                            b := 50

while S <= b - a or E <= abs(f(a)) and E <= abs(f(b)) do
    p := (a + b)/2;
    if f(p) = 0 then
        break;
    elif f(a)*f(p) < 0 then
        b := p;
    else
        a := c;
    end if;
end do;
                            p := 30

Error, cannot determine if this expression is true or false: (950*ln(20/17)-598.10)*(950*ln(20/11)-794.30) < 0



1 Answers

A conditional check inside,

if ... then

tests using Maple's evalb. But your initial values for a and b are exact, not floats, so evalb is not the right check. That's what that error message means, that it cannot tell whether that inequality holds using only the default evalb checker (boolean test).

The most straightforward thing to do is to wrap the expressions on each side of your inequalities (used in the conditional checks) with a call to evalf. Then both sides of the inquality check would become floats, and the regular evalb test would work. That also handles the case that the user-supplied function f itself produces real but exact values instead of just floats.

An alternative could be to utilize is around the inequality in the conditional, but if you don't force any evalf then each iteration of p might become an enormous exact mess. That's not what you want.

You've made at least one other mistake: You have a := c where you mean a := p.

I get an rough idea of what kind of stopping criteria you want, but there are problems with it. Because Maple is using Digits=10 working precision by default then b-a will not become smaller than E. And the check on f(p)=0 is exact and inadequate; it'll likely never become exactly zero. Better might be, say, abs(f(p)) < E or whatever tolerance you decide.

But you really do need stronger stopping criteria. As written it will run away forever. Since the fixed working precision means that you can only subdivide successfully down to a particular floating-point fineness dictated by Digits, you could try something like one of the following.

restart;
E := 1e-8:
f := x -> 950*ln(200000/(200000 - 3000*x)) + (-1)*9.81*x - 500:
a := 10:
b := 50:
while evalf(f(a)*f(b)) < 0 and
      evalf(abs((b-a)/b)) >= 6*10^(-Digits) do
    p := evalf((a + b)/2);
    if abs(f(p)) < E then
        break;
    elif evalf(f(a)*f(p)) < 0 then
        b := p;
    else
        a := p;
    end if;
end do;
a, f(a), p, f(p), b, f(b);

Or, something like, say,

restart;
E := 1e-8:
f := x -> 950*ln(200000/(200000 - 3000*x)) + (-1)*9.81*x - 500:
a := 10:
b := 50:
p := evalf((a + b)/2):
while evalf(f(a)*f(b)) < 0
      and (p-a) <> 0.0 and (b-p) <> 0.0 do
    if abs(f(p)) < E then
        break;
    elif evalf(f(a)*f(p)) < 0 then
        b := p;
    else
        a := p;
    end if;
    p := evalf((a + b)/2);
end do;
a, f(a), p, f(p), b, f(b);

Fell free to put back some of your own stopping crteria, but be careful about brackets. And your need at some guard against p = evalf((a+b)/2) computing at the same as either a or b -- see my two alternates to prevent runaway looping due to this floating-point liability.

Better still, add a loop counter and stop at some maximal allowed number of iterations. You can make it large-ish. Loops that run away forever because of mistakes in the code and lack of defensive programming are a bit ugly.

Also, you didn't actually test that the original values for f(a) and f(b) have opposite signs, prior to entering the loop. You could test that up front, or add it as another stopping criteria (as I've done).

Related