Roots of a cubic polynomial: Difference between revisions
(→{{header|Wren}}: Added note re reversal of coefficients, added all Raku examples and quadratic solved using formula rather than Newton which deals with complex roots.) |
(→{{header|Phix}}: added translation of raku) |
||
Line 23: | Line 23: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Wren}} |
=== {{trans|Wren}} === |
||
Note the seven functions matrix_trace() ... root_poly() are straightforward translations of code from matrix.wren and math.wren, and are not tested beyond matching the output here. There is also some rather suspicious use of reverse() in eigenvalues(), which kinda brute-forces it to work. |
Note the seven functions matrix_trace() ... root_poly() are straightforward translations of code from matrix.wren and math.wren, and are not tested beyond matching the output here. There is also some rather suspicious use of reverse() in eigenvalues(), which kinda brute-forces it to work. |
||
<syntaxhighlight lang="phix"> |
<syntaxhighlight lang="phix"> |
||
Line 227: | Line 227: | ||
errors: {0,0,0} |
errors: {0,0,0} |
||
</pre> |
</pre> |
||
=== {{trans|Raku}} === |
|||
Note the "wrentest" gives eigenvalues of {6,3,-5}... |
|||
<syntaxhighlight lang="phix"> |
|||
with javascript_semantics |
|||
function polynomial(sequence t) |
|||
atom a = 1, // create characteristic polynomial |
|||
b = -(t[1,1] + t[2,2] + t[3,3]), // = -trace |
|||
c = ( t[1,1]*t[2,2] + t[2,2]*t[3,3] + t[3,3]*t[1,1] ) |
|||
-(t[2,3]*t[3,2] + t[3,1]*t[1,3] + t[1,2]*t[2,1] ), |
|||
d = - t[1,1] * t[2,2] * t[3,3] |
|||
- t[1,2] * t[2,3] * t[3,1] |
|||
- t[1,3] * t[2,1] * t[3,2] |
|||
+ t[1,1] * t[3,2] * t[2,3] |
|||
+ t[1,2] * t[2,1] * t[3,3] |
|||
+ t[1,3] * t[2,2] * t[3,1]; // = -determinant |
|||
return {a, b, c, d} |
|||
end function |
|||
function reduction(sequence poly) |
|||
atom {a, b, c, d} = poly, |
|||
b3 = power(b,3), |
|||
b2 = power(b,2), |
|||
a3 = power(a,3), |
|||
a2 = power(a,2), |
|||
delta = 18*a*b*c*d - 4*b3*d + b2*power(c,2) - 4*a*power(c,3) - 27*a2*power(d,2), |
|||
p = (3*a*c - b2) / (3*a2), |
|||
q = (2*b3 - 9*a*b*c + 27*a2*d) / (27*a3), |
|||
d0 = b2 - 3*a*c, |
|||
d1 = 2*b3 - 9*a*b*c + 27*a2*d; |
|||
return { delta, p, q, d0, d1 } |
|||
end function |
|||
include complex.e |
|||
function horner(object x, atom a, b, c, d ) // cubic polynomial using horner's method |
|||
if complex(x) then |
|||
complexn r = 0 |
|||
for f in {a,b,c,d} do |
|||
r = complex_add(complex_mul(r,x),f) |
|||
end for |
|||
return r |
|||
end if |
|||
return ((a*x + b)*x + c)*x + d |
|||
end function |
|||
function solution(atom a, b, /*c*/, /*d*/, delta, p, q, d0, d1) |
|||
sequence x |
|||
if abs(delta)<1e-12 then // multiple real roots |
|||
if abs(p)<1e-12 then // triple equal real roots |
|||
x = {0,0,0} |
|||
else // double real root |
|||
atom x23 = -3*q/(2*p) |
|||
x = {3*q/p, x23, x23} |
|||
end if |
|||
elsif delta > 0 then // three distinct real roots |
|||
x = {} |
|||
for i=1 to 3 do |
|||
x &= 2*sqrt(-p/3) * cos(arccos(sqrt(-3/p)*3*q/(2*p))/3 - 2*PI*(i-1)/3) |
|||
end for |
|||
else // delta < 0 // one real root and two complex conjugate roots |
|||
complexn g |
|||
if d0=0 and d1<0 then |
|||
g = -power(-d1,1/3) |
|||
elsif d0<0 and d1=0 then |
|||
g = -sqrt(-d0) |
|||
else |
|||
complex s = complex_sqrt(power(d1,2) - 4*power(d0,3)), |
|||
g1 = complex_power(complex_div(complex_add(d1,s),2),1/3), |
|||
g2 = complex_power(complex_div(complex_sub(d1,s),2),1/3) |
|||
g = iff(g1={0,0} ? g2 : g1) |
|||
end if |
|||
complex z = complex_mul(g,complex_div(complex_add(-1,complex_sqrt(-3)),2)), |
|||
x1 = complex_mul(-1/3,complex_add(g,complex_div(d0,g))), |
|||
x2 = complex_mul(-1/3,complex_add(z,complex_div(d0,z))) |
|||
x = {x1,x2,complex_conjugate(x2)} |
|||
for i=1 to 3 do |
|||
x[i] = complex_sub(x[i],(1/3)*b/a) |
|||
end for |
|||
return x |
|||
end if |
|||
for i=1 to 3 do |
|||
x[i] -= (1/3)*b/a |
|||
end for |
|||
return x |
|||
end function |
|||
function spectrum(sequence m) |
|||
sequence poly = polynomial(m) |
|||
atom {a, b, c, d} = poly, |
|||
{delta, p, q, d0, d1} = reduction(poly); |
|||
sequence s = solution(a, b, c, d, delta, p, q, d0, d1 ), |
|||
e = apply(true,horner,{s,a,b,c,d}) |
|||
return {s, e} |
|||
end function |
|||
constant r = 1/sqrt(2), |
|||
m = {{"triple", {{ 1, -1, 0}, { 0, 1, -1}, { 0, 0, 1}}}, |
|||
{"wrentest", {{-2, -4, 2}, {-2, 1, 2}, { 4, 2, 5}}}, |
|||
{"double", {{ 2, 0, 0}, { 0, -1, 0}, { 0, 0, -1}}}, |
|||
{"distinct", {{ 2, 0, 0}, { 0, 3, 4}, { 0, 4, 9}}}, |
|||
{"rotation", {{ 1, 0, 0}, { 0, r, -r}, { 0, r, r}}}} |
|||
for i,pair in m do |
|||
{string desc, sequence t} = pair |
|||
printf(1," %d: %8s matrix: %v\n",{i,desc,t}) |
|||
sequence poly = polynomial(t), |
|||
reduct = reduction(poly), |
|||
{s,e} = spectrum(t) |
|||
printf(1," polynomial: %v\n",{poly}) |
|||
printf(1," reduction: %v\n",{reduct}) |
|||
printf(1," eigenvalues: %v\n",{s}) |
|||
printf(1," errors: %v\n",{e}) |
|||
end for |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
1: triple matrix: {{1,-1,0},{0,1,-1},{0,0,1}} |
|||
polynomial: {1,-3,3,-1} |
|||
reduction: {0,0,0,0,0} |
|||
eigenvalues: {1,1,1} |
|||
errors: {0,0,0} |
|||
2: wrentest matrix: {{-2,-4,2},{-2,1,2},{4,2,5}} |
|||
polynomial: {1,-4,-27,90} |
|||
reduction: {69696,-32.33333333,49.25925926,97,1330} |
|||
eigenvalues: {6,3.0,-5.0} |
|||
errors: {0,-4.263256414e-14,-1.705302566e-13} |
|||
3: double matrix: {{2,0,0},{0,-1,0},{0,0,-1}} |
|||
polynomial: {1,0,-3,-2} |
|||
reduction: {0,-3,-2,9,-54} |
|||
eigenvalues: {2,-1,-1} |
|||
errors: {0,0,0} |
|||
4: distinct matrix: {{2,0,0},{0,3,4},{0,4,9}} |
|||
polynomial: {1,-14,35,-22} |
|||
reduction: {8100,-30.33333333,-61.92592593,91,-1672} |
|||
eigenvalues: {11,2.0,1.0} |
|||
errors: {0,1.065814103e-14,-1.77635684e-14} |
|||
5: rotation matrix: {{1,0,0},{0,0.7071067812,-0.7071067812},{0,0.7071067812,0.7071067812}} |
|||
polynomial: {1,-2.414213562,2.414213562,-1.0} |
|||
reduction: {-0.686291501,0.4714045208,-0.0994922778,-1.414213562,-2.686291501} |
|||
eigenvalues: {{1.0,0},{0.7071067812,-0.7071067812},{0.7071067812,0.7071067812}} |
|||
errors: {{-2.22044605e-16,0},{-3.330669073e-16,-3.885780587e-16},{-3.330669073e-16,3.885780587e-16}} |
|||
</pre> |
|||
Complex results on the last entry are just shown as {real,imag} pairs. |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
Revision as of 07:36, 28 November 2023
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]
Phix
Note the seven functions matrix_trace() ... root_poly() are straightforward translations of code from matrix.wren and math.wren, and are not tested beyond matching the output here. There is also some rather suspicious use of reverse() in eigenvalues(), which kinda brute-forces it to work.
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 den = {-roots[1],1},
num = poly_div(reverse(cp),den)
// find second root
roots[2] = root_poly(num)
errs[2] = eval_poly(cp, roots[2])
// divide out to get linear
den = {-roots[2], 1}
num = reverse(poly_div(reverse(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}
Note the "wrentest" gives eigenvalues of {6,3,-5}...
with javascript_semantics
function polynomial(sequence t)
atom a = 1, // create characteristic polynomial
b = -(t[1,1] + t[2,2] + t[3,3]), // = -trace
c = ( t[1,1]*t[2,2] + t[2,2]*t[3,3] + t[3,3]*t[1,1] )
-(t[2,3]*t[3,2] + t[3,1]*t[1,3] + t[1,2]*t[2,1] ),
d = - t[1,1] * t[2,2] * t[3,3]
- t[1,2] * t[2,3] * t[3,1]
- t[1,3] * t[2,1] * t[3,2]
+ t[1,1] * t[3,2] * t[2,3]
+ t[1,2] * t[2,1] * t[3,3]
+ t[1,3] * t[2,2] * t[3,1]; // = -determinant
return {a, b, c, d}
end function
function reduction(sequence poly)
atom {a, b, c, d} = poly,
b3 = power(b,3),
b2 = power(b,2),
a3 = power(a,3),
a2 = power(a,2),
delta = 18*a*b*c*d - 4*b3*d + b2*power(c,2) - 4*a*power(c,3) - 27*a2*power(d,2),
p = (3*a*c - b2) / (3*a2),
q = (2*b3 - 9*a*b*c + 27*a2*d) / (27*a3),
d0 = b2 - 3*a*c,
d1 = 2*b3 - 9*a*b*c + 27*a2*d;
return { delta, p, q, d0, d1 }
end function
include complex.e
function horner(object x, atom a, b, c, d ) // cubic polynomial using horner's method
if complex(x) then
complexn r = 0
for f in {a,b,c,d} do
r = complex_add(complex_mul(r,x),f)
end for
return r
end if
return ((a*x + b)*x + c)*x + d
end function
function solution(atom a, b, /*c*/, /*d*/, delta, p, q, d0, d1)
sequence x
if abs(delta)<1e-12 then // multiple real roots
if abs(p)<1e-12 then // triple equal real roots
x = {0,0,0}
else // double real root
atom x23 = -3*q/(2*p)
x = {3*q/p, x23, x23}
end if
elsif delta > 0 then // three distinct real roots
x = {}
for i=1 to 3 do
x &= 2*sqrt(-p/3) * cos(arccos(sqrt(-3/p)*3*q/(2*p))/3 - 2*PI*(i-1)/3)
end for
else // delta < 0 // one real root and two complex conjugate roots
complexn g
if d0=0 and d1<0 then
g = -power(-d1,1/3)
elsif d0<0 and d1=0 then
g = -sqrt(-d0)
else
complex s = complex_sqrt(power(d1,2) - 4*power(d0,3)),
g1 = complex_power(complex_div(complex_add(d1,s),2),1/3),
g2 = complex_power(complex_div(complex_sub(d1,s),2),1/3)
g = iff(g1={0,0} ? g2 : g1)
end if
complex z = complex_mul(g,complex_div(complex_add(-1,complex_sqrt(-3)),2)),
x1 = complex_mul(-1/3,complex_add(g,complex_div(d0,g))),
x2 = complex_mul(-1/3,complex_add(z,complex_div(d0,z)))
x = {x1,x2,complex_conjugate(x2)}
for i=1 to 3 do
x[i] = complex_sub(x[i],(1/3)*b/a)
end for
return x
end if
for i=1 to 3 do
x[i] -= (1/3)*b/a
end for
return x
end function
function spectrum(sequence m)
sequence poly = polynomial(m)
atom {a, b, c, d} = poly,
{delta, p, q, d0, d1} = reduction(poly);
sequence s = solution(a, b, c, d, delta, p, q, d0, d1 ),
e = apply(true,horner,{s,a,b,c,d})
return {s, e}
end function
constant r = 1/sqrt(2),
m = {{"triple", {{ 1, -1, 0}, { 0, 1, -1}, { 0, 0, 1}}},
{"wrentest", {{-2, -4, 2}, {-2, 1, 2}, { 4, 2, 5}}},
{"double", {{ 2, 0, 0}, { 0, -1, 0}, { 0, 0, -1}}},
{"distinct", {{ 2, 0, 0}, { 0, 3, 4}, { 0, 4, 9}}},
{"rotation", {{ 1, 0, 0}, { 0, r, -r}, { 0, r, r}}}}
for i,pair in m do
{string desc, sequence t} = pair
printf(1," %d: %8s matrix: %v\n",{i,desc,t})
sequence poly = polynomial(t),
reduct = reduction(poly),
{s,e} = spectrum(t)
printf(1," polynomial: %v\n",{poly})
printf(1," reduction: %v\n",{reduct})
printf(1," eigenvalues: %v\n",{s})
printf(1," errors: %v\n",{e})
end for
- Output:
1: triple matrix: {{1,-1,0},{0,1,-1},{0,0,1}} polynomial: {1,-3,3,-1} reduction: {0,0,0,0,0} eigenvalues: {1,1,1} errors: {0,0,0} 2: wrentest matrix: {{-2,-4,2},{-2,1,2},{4,2,5}} polynomial: {1,-4,-27,90} reduction: {69696,-32.33333333,49.25925926,97,1330} eigenvalues: {6,3.0,-5.0} errors: {0,-4.263256414e-14,-1.705302566e-13} 3: double matrix: {{2,0,0},{0,-1,0},{0,0,-1}} polynomial: {1,0,-3,-2} reduction: {0,-3,-2,9,-54} eigenvalues: {2,-1,-1} errors: {0,0,0} 4: distinct matrix: {{2,0,0},{0,3,4},{0,4,9}} polynomial: {1,-14,35,-22} reduction: {8100,-30.33333333,-61.92592593,91,-1672} eigenvalues: {11,2.0,1.0} errors: {0,1.065814103e-14,-1.77635684e-14} 5: rotation matrix: {{1,0,0},{0,0.7071067812,-0.7071067812},{0,0.7071067812,0.7071067812}} polynomial: {1,-2.414213562,2.414213562,-1.0} reduction: {-0.686291501,0.4714045208,-0.0994922778,-1.414213562,-2.686291501} eigenvalues: {{1.0,0},{0.7071067812,-0.7071067812},{0.7071067812,0.7071067812}} errors: {{-2.22044605e-16,0},{-3.330669073e-16,-3.885780587e-16},{-3.330669073e-16,3.885780587e-16}}
Complex results on the last entry are just shown as {real,imag} pairs.
Raku
(as recovered from pdf uploads)
my \r = 1/sqrt(2);
my Pair @m = [[
triple => [[ 1, -1, 0], [ 0, 1, -1], [ 0, 0, 1] ],
double => [[ 2, 0, 0], [ 0, -1, 0], [ 0, 0, -1] ],
distinct => [[ 2, 0, 0], [ 0, 3, 4], [ 0, 4, 9] ],
rotation => [[ 1, 0, 0], [ 0, r, -r], [ 0, r, r] ],
]];
for ^@m {
my Pair $pair = @m[$_]; print "\n {1+$_}: {$pair.key} ";
my @t = $pair.value.[]; say 'matrix: ', @t.raku;
my @poly = polynomial( @t ); say 'polynomial: ', @poly.raku;
my @reduction = reduction( |@poly ); say 'reduction : ', @reduction.raku;
my ( $s, $e ) = spectrum( @t ); say 'eigenvalues: ', @$s.raku;
say 'errors: ', @$e.raku;
}
sub horner( \x, \a, \b, \c, \d ) { # cubic polynomial using horner's method
((a * x + b) * x + c) * x + d
}
sub polynomial( @t ) {
my \a = 1; # create characteristic polynomial
my \b = -(@t[0;0] + @t[1;1] + @t[2;2]); # = -trace
my \c = ( @t[0;0]*@t[1;1] + @t[1;1]*@t[2;2] + @t[2;2]*@t[0;0] )
-(@t[1;2]*@t[2;1] + @t[2;0]*@t[0;2] + @t[0;1]*@t[1;0] );
my \d = - @t[0;0] * @t[1;1] * @t[2;2]
- @t[0;1] * @t[1;2] * @t[2;0]
- @t[0;2] * @t[1;0] * @t[2;1]
+ @t[0;0] * @t[2;1] * @t[1;2]
+ @t[0;1] * @t[1;0] * @t[2;2]
+ @t[0;2] * @t[1;1] * @t[2;0]; # = -determinant
return [ a, b, c, d]>>.narrow;
}
sub reduction( \a, \b, \c, \d ) {
my \delta = 18 * a * b * c * d - 4 * b **3 * d + b**2 * c**2 - 4 * a * c**3 - 27 * a**2 * d**2;
my \p = (3 * a * c - b * b) / (3 * a * a);
my \q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / (27 * a ** 3);
my \d0 = b*b - 3 * a * c;
my \d1 = 2*b**3 - 9*a*b*c + 27 * a**2 * d;
return [ delta, p, q, d0, d1 ];
}
sub solution( \a, \b, \c, \d, \delta, \p, \q, \d0, \d1 ) {
my @x;
if abs(delta) =~= 0 { # say " multiple real roots ", p.raku ;
if abs(p) =~= 0 { # say " triple equal real roots: ";;
@x[$_] = 0 for ^3;
}
else { # say ' double real root:';
@x[0] = 3 * q / p;
@x[1] = -3 * q /(2 * p);
@x[2] = @x[1];
}
}
elsif delta > 0 { # say " three distinct real roots:";
@x[$_] = 2*sqrt(-p/3) * cos( acos( sqrt( -3 / p ) * 3 * q /( 2 * p ) )/3 - 2*pi * $_/ 3 ) for ^3;
}
else { # delta < 0; say " one real root and two complex conjugate roots:";
my $g = do {
if d0 == 0 and d1 < 0 {
-(-d1)**⅓;
}
elsif d0 < 0 and d1 == 0 {
-sqrt(-d0);
}
else {
my \s = Complex( d1**2 - 4 * d0**3 ).sqrt;
my \g1 = (( d1 + s )/2)**⅓;
my \g2 = (( d1 - s )/2)**⅓;
g1 == 0 ?? g2 !! g1
}
}
my Complex $z = $g * ( -1 + Complex(-3).sqrt )/2;
@x[0] = -⅓ * ( $g + d0 / $g );
@x[1] = -⅓ * ( $z + d0 / $z );
@x[2] = @x[1].conj;
}
@x[$_] -= ⅓ * b / a for ^3;
return @x;
}
sub spectrum( @mat ) {
my ( \a, \b, \c, \d ) = polynomial( @mat );
my (\delta, \p, \q, \d0, \d1) = reduction( a, b, c, d );
my @s = solution( a, b, c, d, delta, p, q, d0, d1 )>>.narrow;
my @e = @s.map( { horner( $_, a, b, c, d ) })>>.narrow;
return ( @s, @e );
}
- Output:
1: triple matrix: [[1, -1, 0], [0, 1, -1], [0, 0, 1]] polynomial: [1, -3, 3, -1] reduction : [0, 0.0, 0.0, 0, 0] eigenvalues: [1, 1, 1] errors: [0, 0, 0] 2: double matrix: [[2, 0, 0], [0, -1, 0], [0, 0, -1]] polynomial: [1, 0, -3, -2] reduction : [0, -3.0, -2.0, 9, -54] eigenvalues: [2, -1, -1] errors: [0, 0, 0] 3: distinct matrix: [[2, 0, 0], [0, 3, 4], [0, 4, 9]] polynomial: [1, -14, 35, -22] reduction : [8100, <-91/3>, <-1672/27>, 91, -1672] eigenvalues: [11, 2, 0.9999999999999991e0] errors: [0, 0, -1.0658141036401503e-14] 4: rotation matrix: [[1, 0, 0], [0, 0.7071067811865475e0, -0.7071067811865475e0], [0, 0.7071067811865475e0, 0.7071067811865475e0]] polynomial: [1, -2.414213562373095e0, 2.414213562373095e0, -0.9999999999999998e0] reduction : [-0.6862915010152442e0, 0.4714045207910316e0, -0.0994922778153789e0, -1.414213562373095e0, -2.68629150101523e0] eigenvalues: [0.9999999999999992e0, <0.7071067811865478-0.7071067811865472i>, <0.7071067811865478+0.7071067811865472i>] errors: [0, <-3.3306690738754696e-16-3.885780586188048e-16i>, <-3.3306690738754696e-16+3.885780586188048e-16i>]
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 cubic polynomial after the first root is found and then code from the Roots_of_a_quadratic_function#Wren task to solve the resulting quadratic polynomial which may have complex roots. Note that the former code assumes that polynomials are presented starting from the lowest order term and all other functions used here from the highest order term. It is therefore necessary to reverse the order of terms when switching between the two.
import "./matrix" for Matrix
import "./math" for Math
import "./complex" for Complex
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])
}
}
// assumes real coefficients, highest to lowest order
var quadratic = Fn.new { |a, b, c|
var d = b*b - 4*a*c
if (d == 0) {
// single real root
var root = -b/(a*2)
return [root, root]
}
if (d > 0) {
// two real roots
var sr = d.sqrt
d = (b < 0) ? sr - b : -sr - b
return [d/(2*a), 2*c/d]
}
// two complex roots
var den = 1 / (2*a)
var t1 = Complex.new(-b*den, 0)
var t2 = Complex.new(0, (-d).sqrt * den)
return [t1+t2, t1-t2]
}
// same assumption as quadratic but 'x' can be real or complex.
var evalPoly = Fn.new { |coefs, x|
var c = coefs.count
if (c == 1) return coefs[0]
var sum = coefs[0]
for (i in 1...c) sum = x * sum + coefs[i]
return sum
}
var eigenvalues = Fn.new { |m|
var roots = []
var errs = []
// 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.add(Math.rootPoly(cp))
errs.add(evalPoly.call(cp, roots[0]))
// divide out to get quadratic
var num = Polynom.new(cp[-1..0])
var den = Polynom.new([-roots[0], 1])
var qdr = (num/den)[0]
// find second and third roots
var coefs = qdr.factors[-1..0]
roots = roots + quadratic.call(coefs[0], coefs[1], coefs[2])
errs.add(evalPoly.call(cp, roots[1]))
errs.add(evalPoly.call(cp, roots[2]))
return [cp, roots, errs]
}
var v = 0.70710678118655
var arrs = [
[
[ 1, -1, 0],
[ 0, 1, 1],
[ 0, 0, 1]
],
[
[-2, -4, 2],
[-2, 1, 2],
[ 4, 2, 5]
],
[
[ 1, -1, 0],
[ 0, 1, -1],
[ 0, 0, 1]
],
[
[ 2, 0, 0],
[ 0, -1, 0],
[ 0, 0, -1]
],
[
[ 2, 0, 0],
[ 0, 3, 4],
[ 0, 4, 9]
],
[
[ 1, 0, 0],
[ 0, v, -v],
[ 0, v, v]
]
]
for (i in 0...arrs.count) {
var m = Matrix.new(arrs[i])
System.print("For matrix:")
if (i < arrs.count-1) Fmt.mprint(m, 2, 0) else Fmt.mprint(m, 17, 14)
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, 6, -5] and the corresponding errors are: [0, 0, 0] 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 0 0| | 0 -1 0| | 0 0 -1| whose characteristic polynomial is: x³ - 3x - 2 Its eigenvalues are: [-1, 2, -1] and the corresponding errors are: [0, 0, 0] For matrix: | 2 0 0| | 0 3 4| | 0 4 9| whose characteristic polynomial is: x³ - 14x² + 35x - 22 Its eigenvalues are: [1, 11, 2] and the corresponding errors are: [0, 0, 0] For matrix: | 1.00000000000000 0.00000000000000 0.00000000000000| | 0.00000000000000 0.70710678118655 -0.70710678118655| | 0.00000000000000 0.70710678118655 0.70710678118655| whose characteristic polynomial is: x³ - 2.4142135623731x² + 2.4142135623731x - 1 Its eigenvalues are: [1, 0.70710678118655 + 0.70710678118655i, 0.70710678118655 - 0.70710678118655i] and the corresponding errors are: [0, 0 + 0i, 0 + 0i]