1

Question: Given a cubic equation where the coefficients are real but can take any extreme conditions (e.g very large or very small number), write a program in Matlab that finds all the roots of this equation. You can't use the built-in functions roots and fzero

Here's my approach: By the Fundamental Theorem of Algebra, f(x) must have at least one real root. If I can find this root, I can factor it out and find the other two roots by using the quadratic formula. But I stuck badly on finding the first real root. I'm thinking about the bisection method, but how can I even find the interval that contains the first root. Can anyone help me at this step?

Thanks!

Huy Truong
  • 327
  • Well, since you are programming there are cubic root formulas: https://en.wikipedia.org/wiki/Cubic_function#Roots_of_a_cubic_function – fleablood Oct 08 '15 at 01:37
  • Can you use eig? – Robert Israel Oct 08 '15 at 04:07
  • @RobertIsrael Cute trick :) – Ian Oct 08 '15 at 04:25
  • @RobertIsrael I suppose I can but that is too advanced for this type of question. However, could you give me a hint on how to use "eig"? I don't see any square matrix here to apply it. – Huy Truong Oct 08 '15 at 15:03
  • @fleablood I think just using a straightforward formula will help. There will be a lot of round off errors. – Huy Truong Oct 08 '15 at 15:05
  • @HuyTruong Given a monic polynomial (i.e. leading coefficient $1$), construct its companion matrix (https://en.wikipedia.org/wiki/Companion_matrix), which has that polynomial as its characteristic polynomial. The eigenvalues of the matrix are then the roots of the polynomial. – Robert Israel Oct 08 '15 at 15:47
  • This "trick" actually makes good sense if you have efficient numerical methods for finding eigenvalues, which Matlab has. It is (or was) what Matlab actually uses for finding roots of polynomials. – Robert Israel Oct 08 '15 at 15:49
  • Yes, Matlab documentation says this is what roots does. – Robert Israel Oct 08 '15 at 17:01
  • @RobertIsrael I've just looked up how to construct a companion matrix. If this is allowed, it will take me about 20 lines of code to finish this assignment... – Huy Truong Oct 08 '15 at 22:19

3 Answers3

1

Let's assume the coeff of $x^3$ is positive; the negative case is similar.

The derivative is a quadratic, whose roots $r_1 < r_2$ you can find. (If it has no roots, see below).

To the right of $x = r_2$, the function is increasing, so if $f(r^2) > 0$, then there's no root there, and if $f(r_2) < 0$, then there's at most one root there. A similar argument applies to $f(r_1)$. So we have three cases:

  1. $f(r_1) < 0, f(r_2) > 0$: there's a root between them.

  2. $f(r_1)$ and $f(r_2)$ are both positive: there's a root to the left of $r_1$.

  3. $f(r_1)$ and $f(r_2)$ are both negative: there's a root to the right of $r_2$.

Here's what you can do in case 3: look at $f(r_2 + 1)$, $f(r_2 + 2)$, $f(r_2 + 4)$, ... $f(r_2 + 2^n)$. Eventually one of these will be positive. Then between $r^2$ and this point, there's a root, which you can find by the intermediate value theorem.

A similar thing works for case 2. Case 1 is easy.

And what if the roots of the derivative are both complex? Then the slope is everywhere positive. Look at the points $x =\pm 1$, $x = \pm 2$, ..., $x = \pm 2^n$, increasing $n$ until the signs differ. Now use bisection to find the root.

John Hughes
  • 93,729
  • 1
    A useful trick to save a lot of difficulty: $|x^3+ax^2+bx+c| \geq |x|^3-|a||x|^2-|b||x|-|c|$. Thus any roots must occur for $|x| \leq |a|+|b|+|c|$. (This holds whether the roots are real or complex.) Thus in fact we can do bisection with no brute force at all: just bisect on the interval $[-|a|-|b|-|c|,|a|+|b|+|c|]$. – Ian Oct 08 '15 at 02:08
  • Really nice point. I sure wish I'd thought of that! – John Hughes Oct 08 '15 at 02:11
  • Ah, sorry, that should be maximum with $1$ (otherwise replacing $|x|$ and $1$ with $|x|^2$ achieves the wrong thing. I'm actually going to write an answer to this effect for clarity. – Ian Oct 08 '15 at 02:14
  • By the way, the reason I remember this is because one can write what I just said and apply Rouche's theorem to get a nice proof of the Fundamental Theorem of Algebra. :) – Ian Oct 08 '15 at 02:36
  • Even better. And I thought that I'd seen just about every possible application of Rouche's theorem -- clearly not. :) – John Hughes Oct 08 '15 at 11:45
1

This may help: If $x^3 +A x^2+B x+C=0$ then $$|x|<1+\max (|A|,|B|,|C|).$$ And if $C\ne 0$ then $$|x|>(1+\max (|B/C|,|A/C|,|1/C|))^{-1}.$$ Sometimes this can be sharpened by letting $x=k y$ with $k>0$ because we have $$|x|<k(1+\max (|A/k|,|B/k^2|,|C/k^3|)=k+\max (|A|,|B/k|,|C/k^2|),$$ with a corresponding lower bound for |x| in terms of $A,B,C,k$ from the 2nd inequality. These inequalities are verifiable by entirely elementary algebra.

DanielWainfleet
  • 57,985
  • Also in the case $C=0$ the situation is rather simple. – Ian Oct 08 '15 at 02:56
  • In general if $ p(x)=\sum_{j=0}^{j=n} A_j x^j$ with $ A_n=1$ and $n>0$ we have $ p(x)=0 \to |x|<1+\max (|A_j| : j<n)$. If I wanted to be really pedantic I would add "with the convention that $0^0=1$" to define $A_0 x^0$ when $x=0$. – DanielWainfleet Oct 08 '15 at 04:07
0

Normalize so that the leading coefficient is $1$. By the triangle inequality:

$$|x^3+ax^2+bx+c| \geq |x|^3-|a||x|^2-|b||x|+|c|.$$

Now let $R=\max \{ 1,|a|+|b|+|c| \}$. Then the inequality above implies that all roots must lie in $[-R,R]$. Thus we can find one of the real roots (or the only one) to within an accuracy of $\varepsilon$ with essentially $\log_2(R/\varepsilon)$ bisection steps with the initial interval $[-R,R]$.

Now this might be too slow, if $R$ is big or $\varepsilon$ is small. Indeed, this simple method is not "secure", in the sense that I can cook up a problem which is within the domain of your problem statement and which will make this method take an arbitrarily long time to run. Yet somehow these problems with huge $R$ should actually be very easy.

And indeed, better methods exist which can handle large $R$. As usual it is helpful to know something analytical about the problem in advance. In this case we can find the extrema explicitly using the quadratic formula. If none exist, then the cubic is increasing, which means we can apply Newton's method starting anywhere we want. It is impossible for only one extremum to exist in a cubic (why?) If two extrema $r_1 < r_2$ exist, then you can compute the value of the cubic at the two extrema. If the signs differ, then the root is in between them, so you can run Newton's method starting at $\frac{r_1+r_2}{2}$. If they are both negative, then the root is to the right of them, so you can run Newton's method starting at $\frac{r_2+R}{2}$. If they are both positive, then the root is to the left of them, so you can run Newton's method starting at $\frac{r_1-R}{2}$.

As an example, with $a=-169,b=45,c=-29$, we have two extrema, one around $0$ and one at about $110$. Both are negative and $R$ is about $250$, so we start Newton's method at about $180$ and get a root of about $169$ in four steps. Bisection would have required 40-50 steps.

Ian
  • 101,645
  • I don't get why that inequality implies the roots are in $[-R,R]$. Can you explain that? – bartgol Oct 08 '15 at 03:09
  • 1
    @bartgol For $|x| \geq 1$, $|x^3+ax^2+bx+c| \geq |x|^3-|a||x|^2-|b||x|-|c| \geq |x|^3-|x|^2(|a|+|b|+|c|)$. Factoring, the previous expression is $|x|^2(|x|-|a|-|b|-|c|)$. If $|x|>|a|+|b|+|c|$ then this last expression is strictly positive (assuming of course that we don't have $a=b=c=0$ in which case the situation is extremely simple). – Ian Oct 08 '15 at 03:28
  • Ok, I see it now, thanks. Btw, I think you should have $-|c|$ on the right of the first inequality, no? – bartgol Oct 08 '15 at 03:34
  • @bartgol Yes, I fixed that already. :) – Ian Oct 08 '15 at 03:35
  • @Ian I'm not sure how I can implement your method to find root of equations like $$x^3$$. Such equations have only one root and Newton's method will fail to find the root. Am I correct? – Huy Truong Oct 17 '15 at 04:48
  • @HuyTruong On the contrary, in the monotone case Newton's method will always work. The convergence will be slow because of the multiple root, but it will still converge. – Ian Oct 17 '15 at 13:24
  • @Ian I'm not sure if we are making things extra complicated, but I wrote a very simple code for newton's method (below) and it always successfully gives me the 1st root, no matter what p0 I choose. Can you explain why?

    function x = newton(C, p0, TOL) f = @(x) C(1)*x^3 + C(2)*x^2 + C(3)*x + C(4); df = @(x) 3*C(1)*x^2 + 2*C(2)*x + C(3); while 1 p = p0 - f(p0)/df(p0); if abs(p - p0) < TOL x = p; break end p0 = p; end end

     – Huy Truong Oct 20 '15 at 04:11
  • @Ian Also, for $$x^3$$, it will not work for starting point x = 0. So why does Newton's method work then? – Huy Truong Oct 20 '15 at 04:35
  • @HuyTruong If you start at the root, you need to immediately stop, not try to compute a Newton's method step. So your break condition should actually be at the top of the loop. I'm not sure why your code appears to work independent of initial condition. In the case where you have extrema, I would expect some difficulties if you start in between an extremum and a root, with no other roots on the other side (so that you have to get past the extremum in order to proceed). Such a thing would be like $-x(x+1)^2+0.1$ for initial conditions less than $-1$. – Ian Oct 20 '15 at 12:47
  • @HuyTruong Indeed, if I start at $-3$ with the function I wrote above, MATLAB's fsolve (with default parameters) gives up. – Ian Oct 20 '15 at 12:53
  • @Ian My code have worked perfectly so far for most of the cases I can think of, even for the cubic equation you suggested above. The only difficulty is the case of $$x^3$$. In fact, it gives the answer, ans = 1.0e-09 * 0.1592 + 0.0000i -0.0796 + 0.1379i -0.0796 - 0.1379i when I set tolerance = 1e-10. I can easily modify this in my code by adding a specific case for x^3, but I want to understand why Newton's method doesn't work in this case. – Huy Truong Oct 20 '15 at 13:56
  • That doesn't make any sense...your iteration, when simplified, is $x_{k+1}=\frac{2x_k}{3}$. There should be no complex numbers involved. The error size doesn't surprise me, however, because you are getting subtraction of nearly equal numbers in this situation. You can fix that by simplifying your expression with a common denominator: $x-\frac{ax^3+bx^2+cx+d}{3ax^2+2bx+c}=\frac{3ax^3+2bx^2+cx-ax^3-bx^2-cx-d}{3ax^2+2bx+c}=\frac{2ax^3+bx^2}{3ax^2+2bx+c}$. – Ian Oct 20 '15 at 14:04
  • @HuyTruong Sorry, I made an error: the $-d$ does not cancel, so you have $\frac{2ax^3+bx^2-d}{3ax^2+2bx+c}$. Note that this might not always be better, it depends on the relative size of the numbers. But when $a$ is much larger than everything else, this is definitely the way to go. – Ian Oct 20 '15 at 14:15
  • @Ian Thanks. But yeah, it doesn't makes sense to me either... Here's the way I can explain why there are complex numbers involve: I want to find the first real root, then I factor out my polynomial with deconv, leaving a quadratic polynomial. Then I use quadratic formula to find all the other roots. But since my 1st root is not accurate, the quadratic formula is not accurate, and hence the complex numbers involve. Man, that's really all I can think of, because I've tried with all others cases, and it just works perfectly. – Huy Truong Oct 20 '15 at 15:56
  • @HuyTruong Oh, you're only finding one root and then factoring. Then I think everything makes sense: you obtain a small positive root instead of zero. So when you divide, the quadratic you get is not exactly $x^2$. A small perturbation of $x^2$ can easily have complex roots, such as $x^2+\delta$ for small $\delta>0$. You can reduce this problem by using the roots of the estimated quadratic as initial guesses for a (perhaps complex) Newton iteration. – Ian Oct 20 '15 at 16:01