Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl: the Markov chain saw
 
PerlMonks  

solve cubic equations

by no_slogan (Deacon)
on May 03, 2017 at 05:51 UTC ( [id://1189383]=CUFP: print w/replies, xml ) Need Help??

Everybody knows the quadratic formula, which lets you solve this equation: a x2 + b x + c = 0. Turns out it's not hard to solve when there's also an x3 term. There are either one or three solutions. This algorithm makes me happy.
use constant pi => 3.141592653589793; sub cubic { # solve a cubic equation in the form # x^3 + a x^2 + b x + c = 0 my ($a, $b, $c) = @_; my $q = $a*$a/9 - $b/3; my $r = ($a*$a/27 - $b/6)*$a + $c/2; my $s = $a / -3; my $d = $r*$r - $q*$q*$q; if ($d > 0) { my $t = (sqrt($d) + abs($r)) ** (1/3); my $u = ($t + $q / $t); return $r > 0 ? $s - $u : $s + $u; } else { my $t = atan2(sqrt(-$d), $r) / 3; my $u = 2 * sqrt($q); # $d <= 0 implies $q >= 0 return ( $s - $u * cos($t), $s - $u * cos($t + 2/3*pi), $s - $u * cos($t - 2/3*pi), ); } }

Replies are listed 'Best First'.
Re: solve cubic equations
by vrk (Chaplain) on May 03, 2017 at 08:02 UTC

    It's always fun to implement numerical algorithms, especially when a closed form solution exists. You're missing some roots, though... The fundamental theorem of algebra shows that a polynomial of degree n has exactly n roots, although some (or all) of them can be complex. Your cubic solver finds only the real roots.

    By the way, for finding the roots of polynomials of higher degree, there's an iterative solver on CPAN (Math::Polynomial::Solve). You'll need it if you go higher than degree 5.

      Quite so. If you want the complex roots, use $t multiplied by the three cube roots of unity (1 and -1/2 ± i*sqrt(3)/2) in the first if branch.
Re: solve cubic equations
by no_slogan (Deacon) on May 03, 2017 at 16:33 UTC
    BTW, the some of the solutions returned by the second if branch may be the same. Mathematicians call these "multiple roots." If $u == 0, there is a triple root at $x = $s. If $t == 0, there is a single root at $x = $s - $u and a double root at x = $s + $u/2 (because cos(2/3*pi) == cos(-2/3*pi) == -1/2).
Re: solve cubic equations
by bduggan (Pilgrim) on May 05, 2017 at 19:30 UTC
    I prefer the solution involving radicals.

    What's cool is that you can use square roots to solve quadratics, third roots to solve cubics, fourth roots to solve quartics, but you can't use just radicals from 5 and up.

    Perl 6's built-in roots function comes in handy.

    This equation comes right from the wikipedia entry, and uses the names of the variables there.
    #!/usr/bin/env perl6
    
    sub cubic(\a,\b,\c,\d) {
        my \Δ0 = b²	- 3 × a × c;
        # note: special case when Δ0 == 0
        my \Δ1 = 2 * b³ - 9 × a × b × c + 27 × a² × d;
        my \C = ( ( Δ1 + sqrt( Δ1² - 4 × Δ0³ + 0i) ) / 2 ).roots(3)[0];
        my \ς = 1.roots(3);  # cubic roots of unity
        return [0,1,2].map: -> \k {
            ( -1 / ( 3 × a ) ) × ( b + ς[k] × C + Δ0 / ( C × ς[k] ) )
        }
    }
    
    
    my @vals = cubic(1,10,10,-10);
    
    # test
    use Test;
    plan 3;
    
    my $f = -> \x { x³ + 10 * x² + 10 * x - 10 };
    
    is-approx $f( @vals[0] ), 0, 'first value';
    is-approx $f( @vals[1] ), 0, 'second value';
    is-approx $f( @vals[2] ), 0, 'third value';
    
    
    I had a hard time doing this with a code block on perlmonks, so I made a gist instead.
      my \C = ( ( Δ1 + sqrt( Δ1² - 4 × Δ0³ + 0i) ) / 2 ).roots(3)[0];

      Mathematically, your code is equivalent. (The arctan and cos are hidden inside roots.) Numerically, yours is potentially unstable. The first + in the line above might lose accuracy due to subtracting nearly-equal numbers. See Numerical Recipes, chapter 5.6. Unfortunately, sometimes we have to choose between code that looks nice and code that works.

      The trouble with Perl 6 is that it's full of people who think it's a good idea to have a class called "Cool" which contains a hodgepodge of numerical and string-manipulation methods.

        1. Yes, it's totally possible to optimize for accuracy, and yes, some situations call for that.
        2. Allomorphic typing is cool.
Re: solve cubic equations (Python)
by Anonymous Monk on May 03, 2017 at 15:58 UTC
    This is an area where I actually prefer to use Python instead of Perl, anytime I implement MATHS!
    import math def cubic(a,b,c): q = a*a/9 - b/3 r = (a*a/27 - b/6)*a + c/2 s = a/-3 d = r*r - q*q*q if d > 0: t = (math.sqrt(d) + abs(r)) ** (1/3) u = (t + q/t) return s - u if r > 0 else s + u else: t = math.atan2(math.sqrt(-d), r) / 3 u = 2 * math.sqrt(q) return ( s - u * math.cos(t), s - u * math.cos(t + 2/3*math.pi), s - u * math.cos(t - 2/3*math.pi) ) print( cubic(10,10,-10) )
      Not sure what the advantage of Python is in this case. Your code is basically the same, except that it doesn't have any dollar signs or semicolons. I never liked Python's "x if cond else y" as a replacement for "cond ? x : y", it just doesn't seem right to put the condition in the middle of the expression. But I have to admit that Python's built-in multiple precision ints are a huge advantage for certain kinds of math programming.
        "Your code is basically the same, except that it doesn't have any dollar signs or semicolons."

        Precisely why it looks more like MATHS. :)

          A reply falls below the community's threshold of quality. You may see it by logging in.
      This code does not produce the right answer.
      It produces
      (-9.433363736780528, -9.433363736780528, -9.433363736780528)
      Adding ".0" to all the numbers in the code seems to help.
        Yep that's Python, if the first assignment (declaration) is an integer the variable is typed to be an integer and all operators will only perform integer operations.

        Cheers Rolf
        (addicted to the Perl Programming Language and ☆☆☆☆ :)
        Je suis Charlie!

        Thanks. The point was simply to show how using the right tool for the job leads to a cleaner solution. For all things string related, I would pick Python as a last resort.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: CUFP [id://1189383]
Approved by beech
Front-paged by Arunbear
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others meditating upon the Monastery: (3)
As of 2024-11-10 00:55 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    chatterbot is...






    Results (36 votes). Check out past polls.