Roots of a cubic polynomial

Roots of a cubic polynomial is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Find all eigenvalues of a real 3x3 matrix and estimate their errors.

See Wikipedia Cubic equation. [1] Example.

a*x**3 + b*x**2 + c*x + d = 0

matrix: [[1,-1,0], [0,1,-1],[0,0,1]]

polynomial: [a,b,c,d] =[1, -3,+3, -1]

roots: [1,1,1]

matrix: [[1,-1,0], [0,1,-1],[0,0,1]]

polynomial: [a,b,c,d] =[1, -3,+3, -1]

roots: [1,1,1]

RPL

Works with: HP version 49

« JORDAN 4 ROLL DROP NIP
» 'TASK' STO

[[1 -1 0][0 1 -1][0 0 1]] TASK

Output:

2: 'X^3-3*X^2+3*X-1'
1: [1 1 1]

If the characteristic polynomial is not needed, the EGVL function directly returns the eigenvalues as a vector.

Wren

Library: Wren-matrix

Library: Wren-math

Library: Wren-fmt

The eigenvalues of a 3 x 3 matrix will be the roots of its characteristic polynomial.

We borrow code from the Polynomial_long_division#Wren task to divide out this polynomial after each root is found.

import "./matrix" for Matrix
import "./math" for Math
import "./fmt" for Fmt

class Polynom {
    // assumes factors start from lowest order term
    construct new(factors) {
        _factors = factors.toList
    }

    factors { _factors.toList }

    /(divisor) {
        var curr = canonical().factors
        var right = divisor.canonical().factors
        var result = []
        var base = curr.count - right.count
        while (base >= 0) {
            var res = curr[-1] / right[-1]
            result.add(res)
            curr = curr[0...-1]
            for (i in 0...right.count-1) {
                curr[base + i] = curr[base + i] - res * right[i]
            }
            base = base - 1
        }
        var quot = Polynom.new(result[-1..0])
        var rem = Polynom.new(curr).canonical()
        return [quot, rem]
    }

    canonical() {
        if (_factors[-1] != 0) return this
        var newLen = factors.count
        while (newLen > 0) {
            if (_factors[newLen-1] != 0) return Polynom.new(_factors[0...newLen])
            newLen = newLen - 1
        }
        return Polynom.new(_factors[0..0])
    }
}

var eigenvalues = Fn.new { |m|
    var roots = List.filled(3, 0)
    var errs  = List.filled(3, 0)
    // find the characteristic polynomial
    var cp = List.filled(4, 0)
    cp[0] = 1
    cp[1] = -m.trace
    cp[2] = m.cofactors.trace
    cp[3] = -m.det
    // find first root
    roots[0] = Math.rootPoly(cp)
    errs[0] = Math.evalPoly(cp, roots[0])
    // divide out to get quadratic
    var num = Polynom.new(cp[-1..0])
    var den = Polynom.new([-roots[0], 1])
    num = (num/den)[0]
    // find second root
    roots[1] = Math.rootPoly(num.factors[-1..0])
    errs[1] = Math.evalPoly(cp, roots[1])
    // divide out to get linear
    den = Polynom.new([-roots[1], 1])
    num = (num/den)[0]
    // find third root
    roots[2] = -num.factors[0]
    errs[2] = Math.evalPoly(cp, roots[2])
    return [cp, roots, errs]
}

var mats = [
    Matrix.new ( [
        [ 1, -1,  0],
        [ 0,  1,  1],
        [ 0,  0,  1]
    ]),
    Matrix.new ([
        [-2, -4,  2],
        [-2,  1,  2],
        [ 4,  2,  5]
    ])
]

for (m in mats) {
    System.print("For matrix:")
    Fmt.mprint(m, 2, 0)
    var eigs = eigenvalues.call(m)
    Fmt.print("\nwhose characteristic polynomial is:")
    Fmt.pprint("$n", eigs[0], "", "x")
    Fmt.print("\nIts eigenvalues are: $n", eigs[1])
    Fmt.print("and the corresponding errors are: $n\n", eigs[2])
}

Output:

For matrix:
| 1 -1  0|
| 0  1  1|
| 0  0  1|

whose characteristic polynomial is:
x³ - 3x² + 3x - 1

Its eigenvalues are: [1, 1, 1]
and the corresponding errors are: [0, 0, 0]

For matrix:
|-2 -4  2|
|-2  1  2|
| 4  2  5|

whose characteristic polynomial is:
x³ - 4x² - 27x + 90

Its eigenvalues are: [3, -5, 6]
and the corresponding errors are: [0, 0, 0]