Finding Zeros of cubic without using fzero or roots in Matlab

Question

okay, so I've modifies my code a bit.

function rts = cubicSS(C)
a = C(1);
b = C(2);
c = C(3);
d = C(4);
rts = [90000 ; 90000 ; 90000];
f = @(x) a.*x^3 + b.*x^2 + c.*x + d;
a1=1;
while f(a1) <= 0 && a1 < inf 
    a1=2*a1;
end
b1=1;
while f(-b1)>=0 && b1 < inf
    b1=2*b1;
end
rcount = 0;
tol = 1e-6;
a1=1;
while f(a1) <= 0 && a1 < inf 
    a1=2*a1;
end
b1=1;
while f(-b1)>=0 && b1 < inf
    b1=2*b1;
end

 while 1  
            p = (a1+b1)/2;
            if f(p) == 0
                rcount = rcount + 1;
                rts(rcount) = p;
                break; end
            if p-a1 <tol
               break; end
            if f(a1)*f(p)> 0
                a1=p;
            else
                b1=p;
            end
        end
    end

However, whenever I run this code, I don't get any roots back. right now I'm setting C =[ 1 -3 3 -1] just to test it out. My plan is that after i find one root, I will use deconv function to do long division and then use the quadratic formula on my newly found polynomial .

wanted to know what bugs or errors could possibly be in my code

Answer 1

Since you set your left endpoint as a positive and found its position by taking $f(-b1)$, you start out in the range $[1,2]$ instead of the range $[-1,2]$ like you want. After finding the absolute value of $b1$, make sure it is negative, since that is how your programme works.

Aside from that, you are setting your tolerance based on the distance from $p$ to $a1$, which will eventually terminate, but you never set the value of the root for the "tolerant" version.

I didn't see anything else wrong here.

NOTE, though, that your example cubic will require complex roots.

Answer 2

The code for a1 and b1 in the following is only for cubic polynomial with non-zero leading coefficient!

So, I clean up the code. It was basically right. I've no idea why you get inf for the function values. I change following things

1. I filled rts with nan. Otherwise it was little difficult to read normal size numbers and 90000 the same time.
2. I parenthesis $f$. It is called Horner's method and helps against rounding errors. It is just a rule of thumb and won't always help.
3. b1 is always negative.
4. tol_interval controls, when is the interval from b1 to a1 so small that you consider b1 equals a1.
5. tol_func controls, when is the absolute function value is so small that you consider it zero.

Play with tol_interval and tol_func and you see how ill conditioned finding zeros of polynomial can be. In particular your example, where $1$ is a zero of multiplicity $3$.

Here is the code:

function rts = cubicSS(C)
    a = C(1);
    b = C(2);
    c = C(3);
    d = C(4);
    rts = [nan ; nan ; nan];
    f = @(x) ((a .* x + b) .* x + c) .* x + d;

    s = sign(a);

    a1 = 1; 
    while s*f(a1) <= 0 && a1 < inf
        a1=2*a1;
    end

    b1 = -1;
    while s*f(b1) >= 0 && b1 > -inf
        b1 = 2*b1;
    end

    rcount = 0;
    tol = 1e-6;
    tol_interval = tol*abs(a1 - b1);
    tol_func = tol*max(abs(f(a1)), abs(f(b1)));

    while 1
        p = (a1+b1)/2;
        if abs(f(p)) <= tol_func
            rcount = rcount + 1;
            rts(rcount) = p;
            break; 
        end

        if abs(a1 - b1) <= tol_interval         
            break; 
        end

        if f(a1)*f(p)> 0
            a1=p;
        else
            b1=p;
        end
    end
end

Answer 3

If your polynomial is $f(x) =-x^3-1$, the initial loop on $a1$ will overflow since $f(x) < 0$ for all $x \ge 0$.