Roots of a cubic polynomial
Find all eigenvalues of a real 3x3 matrix and estimate their errors.
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.
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]
Raku
Works with: Raku version 2023.11
Phix
with javascript_semantics
function matrix_trace(sequence m)
atom res = 0
for i=1 to length(m) do
res += m[i,i]
end for
return res
end function
function matrix_minor(sequence m, integer x, y)
integer l = length(m)-1
sequence res = repeat(repeat(0,l),l)
for i=1 to l do
integer ri = i+(i>=x)
for j=1 to l do
integer rj = j+(j>=y)
res[i,j] = m[ri,rj]
end for
end for
return res
end function
function matrix_determinant(sequence m)
// Returns the determinant of the current instance
// provided it's square, using Laplace expansion.
integer l = length(m)
if l=1 then return m[1,1] end if
if l=2 then return m[2,2]*m[1,1]-m[1,2]*m[2,1] end if
integer sgn = 1
atom det = 0
for i=1 to l do
sequence mm = matrix_minor(m,1,i)
det += sgn*m[1,i]*matrix_determinant(mm)
sgn = -sgn
end for
return det
end function
function cofactors(sequence m)
integer l = length(m)
sequence res = repeat(repeat(0,l),l)
for i=1 to l do
for j=1 to l do
sequence mm = matrix_minor(m,i,j)
atom d = matrix_determinant(mm)
res[i,j] = d*power(-1,i+j)
end for
end for
return res
end function
function eval_poly(sequence coefs, atom x)
atom res = coefs[1]
for i=2 to length(coefs) do res = res*x + coefs[i] end for
return res
end function
function diff_poly(sequence coefs)
// Returns the coefficients of a polynomial following differentiation.
// Polynomials are represented as described in 'eval_poly' above(!).
integer c = length(coefs)-1
if c=0 then return {0} end if
sequence deriv = repeat(0,c)
for i=1 to c do deriv[i] = (c-i+1) * coefs[i] end for
return deriv
end function
function root_poly(sequence coefs, atom guess=0.001, tol=1e-15, maxIter=100, mult=1)
integer deg = length(coefs)-1
if deg=0 then return null end if
if deg=1 then return -coefs[2]/coefs[1] end if
if deg=2 and coefs[2]*coefs[2] - 4*coefs[1]*coefs[3] < 0 then return null end if
if eval_poly(coefs, 0)=0 then return 0 end if
sequence deriv = diff_poly(coefs)
atom eps = 0.001,
x0 = guess
for iter=1 to maxIter do
atom den = eval_poly(deriv, x0)
if den=0 then
x0 = iff(x0>=0 ? x0 + eps : x0 - eps)
else
atom num = eval_poly(coefs, x0),
x1 = x0 - num/den * mult
if abs(x1-x0)<=tol then
atom r = round(x1)
if abs(r-x1)<=eps and eval_poly(coefs, r)=0 then return r end if
return x1
end if
x0 = x1
end if
end for
x0 = round(x0)
if eval_poly(coefs, x0)=0 then return x0 end if
return null
end function
function poly_div(sequence num,den)
sequence curr = trim_tail(num,0),
right = trim_tail(den,0),
res = {}
integer base = length(curr)-length(right)
while base>=0 do
atom r = curr[-1] / right[-1]
res &= r
curr = curr[1..-2]
for i=1 to length(right)-1 do
curr[base+i] -= r * right[i]
end for
base -= 1
end while
return res
end function
function poly(sequence si)
-- display helper
string r = ""
for t=length(si) to 1 by -1 do
integer sit = si[t]
if sit!=0 then
if sit=1 and t>1 then
r &= iff(r=""? "":" + ")
elsif sit=-1 and t>1 then
r &= iff(r=""?"-":" - ")
else
if r!="" then
r &= iff(sit<0?" - ":" + ")
sit = abs(sit)
end if
r &= sprintf("%d",sit)
end if
r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
end if
end for
if r="" then r="0" end if
return r
end function
function eigenvalues(sequence m)
// find the characteristic polynomial
sequence cp = {1,
-matrix_trace(m),
matrix_trace(cofactors(m)),
-matrix_determinant(m)},
roots = repeat(0,3),
errs = repeat(0,3)
// find first root
roots[1] = root_poly(cp)
errs[1] = eval_poly(cp, roots[1])
// divide out to get quadratic
// aside: no idea why reverse helps...
sequence num = reverse(cp),
den = {-roots[1],1}
num = reverse(poly_div(num,den))
// find second root
roots[2] = root_poly(reverse(num))
errs[2] = eval_poly(cp, roots[2])
// divide out to get linear
den = {-roots[2], 1}
num = reverse(poly_div(num,den))
// find third root
roots[3] = -num[1]
errs[3] = eval_poly(cp, roots[3])
string pm = ppf(m,{pp_Nest,1,pp_Indent,8,pp_IntFmt,"%2d"})
return {pm, poly(reverse(cp)), roots, errs}
end function
constant tests = {{{ 1, -1, 0},
{ 0, 1, 1},
{ 0, 0, 1}},
{{-2, -4, 2},
{-2, 1, 2},
{ 4, 2, 5}}}
constant fmt = """
Matrix: %s
characteristic polynomial: %s
eigenvalues: %v
errors: %v
"""
for m in tests do
printf(1,fmt,eigenvalues(m))
end for
- Output:
Matrix: {{ 1,-1, 0},
{ 0, 1, 1},
{ 0, 0, 1}}
characteristic polynomial: x^3 - 3x^2 + 3x - 1
eigenvalues: {1,1,1}
errors: {0,0,0}
Matrix: {{-2,-4, 2},
{-2, 1, 2},
{ 4, 2, 5}}
characteristic polynomial: x^3 - 4x^2 - 27x + 90
eigenvalues: {3,-5,6}
errors: {0,0,0}
RPL
« 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
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]