Roots of a cubic polynomial: Difference between revisions
m
→{{header|Phix}}: online tag
(added Phix) |
m (→{{header|Phix}}: online tag) |
||
| (19 intermediate revisions by 6 users not shown) | |||
Line 15:
=={{header|J}}==
Using J's [[j:Vocabulary/pdot|polynomial primitive]] for most of the heavy lifting here. Also using J's [[j:Vocabulary/dot|generalized determinant]] to find the characteristic polynomial (first turning each element of the matrix into a polynomial then finding the determinant of the result). See also [[j:Doc/Articles/Play141|Oh, No, Not Eigenvalues Again!]].
Implementation:
<syntaxhighlight lang=J>
NB. matrix characteristic polynomial
mcp=: {{ ;pdf/ .ppr y,each -=i.#y }}
pdf=: -/@,:L:0 NB. polynomial difference
ppr=: +//.@(*/)L:0 NB. polynomial product
proots=: 1 {:: p. NB. polynomial roots
task=: {{
poly=. mcp y
roots=. proots poly
errors=. poly p."1 0 roots
lines=. (,:'matrix: '),&.|:&": y
lp=. 'coefficients: ',":poly
lr=. 'roots: ',":roots
le=. 'errors: ',":errors
lines,lp,lr,:le
}}</syntaxhighlight>
Task examples:
<syntaxhighlight lang=J>
task 1 _1,0 1 1,:0 0 1
matrix: 1 _1 0
0 1 1
0 0 1
coefficients: 1 _3 3 _1
roots: 1 1 1
errors: 0 0 0
task _2 _4 2,_2 1 2,:4 2 5
matrix: _2 _4 2
_2 1 2
4 2 5
coefficients: _90 27 4 _1
roots: 6 _5 3
errors: 0 0 0
task 1 _1,0 1 _1,:0 0 1
matrix: 1 _1 0
0 1 _1
0 0 1
coefficients: 1 _3 3 _1
roots: 1 1 1
errors: 0 0 0
task 2 0,0 _1,:0 0 _1
matrix: 2 0 0
0 _1 0
0 0 _1
coefficients: 2 3 0 _1
roots: 2 _1 _1
errors: 0 0 0
task 2 0,0 3 4,:0 4 9
matrix: 2 0 0
0 3 4
0 4 9
coefficients: 22 _35 14 _1
roots: 11 2 1
errors: 0 0 0
task 1 0,0,.+.1ad_45 1ad45
matrix: 1 0 0
0 0.707107 _0.707107
0 0.707107 0.707107
coefficients: 1 _2.41421 2.41421 _1
roots: 1 0.707107j0.707107 0.707107j_0.707107
errors: 0 4.44089e_15j_1.77636e_15 4.44089e_15j1.77636e_15</syntaxhighlight>
We could instead use [[j:Vocabulary/LAPACK|lapack]]. The lapack approach would be faster and more stable for very large matrices, but that's not an issue here. Also, the approach we use here better fits the current title of this task.
=={{header|jq}}==
'''rootPoly is adapted from [[#Wren|Wren]]'''
The first few functions are taken from
[[Polynomial_long_division#jq]] but are reproduced here for clarity.
In the following, a method for computing the roots of a polynomial
of arbitrary degree is used, and so no restriction to polynomials of
degree 3 is imposed.
Note that the polynomial of degree n, P = SIGMA (c[i] * x^i), is
represented by the array [ c[0] ... c[n] ], and the estimated error
in a root, r, is computed using the formula (dY / (dP/dX)(r)) except
when the denominator or the expression as a whole is computed as 0,
in which case a very small value (almost indistinguishable from 0.0)
is used.
<syntaxhighlight lang="jq">
# Input: a possibly empty numeric array [c0, ... ]
# Emit the canonical form of the polynomial: SIGMA c[i] * x^i
def canonical:
if length == 0 then [0]
elif .[-1] == 0 then .[:-1]|canonical
else .
end;
# string representation
def poly2s: "Polynomial(\(join(",")))";
# Polynomial division
# Output [ quotient, remainder]
def divrem($divisor):
($divisor|canonical) as $divisor
| { curr: canonical}
| .base = ((.curr|length) - ($divisor|length))
| until( .base < 0;
(.curr[-1] / $divisor[-1]) as $res
| .result += [$res]
| .curr |= .[0:-1]
| reduce range (0;$divisor|length-1) as $i (.;
.curr[.base + $i] += (- $res * $divisor[$i]) )
| .base += -1
)
| [(.result | reverse), (.curr | canonical)] ;
# Evaluate a polynomial at $x
# Input: the array representation of a polynomial
def evalPoly($x):
length as $c
| . as $coeffs
| reduce range(1;length+1) as $i (0; . * $x + $coeffs[$c - $i]) ;
# Differentiate a polynomial
# Input: the array representation of the polynomial to be differentiated
def diffPoly:
(length - 1) as $c
| if $c == 0 then [0]
else . as $coefs
| reduce range(0; $c) as $i ([]; .[$i] = ($i+1) * $coefs[$i+1])
end;
# Emit a stream
# No check is made that the input represents a quadratic polynomial
def quadraticRoots:
. as $coefs
| ($coefs[1]*$coefs[1] - 4*$coefs[0]*$coefs[2]) as $d
| select($d >= 0)
| ($d|sqrt) as $s
| -($s + $coefs[1]), ($s - $coefs[1])
| . / (2 * $coefs[2]);
# Attempt to find a real root of a polynomial using Newton's method from
# an initial guess, a given tolerance, a maximum number of iterations
# and a given multiplicity (usually 1).
# If a root is found, it is returned, otherwise 'null' is returned.
# If the root is near an integer, check if the integer is in fact the root.
# If the polynomical degree is 0, return null.
# Input: [c0, c1 ... ]
def rootPoly($guess; $tol; $maxIter; $mult):
. as $coefs
| (length - 1) as $deg
| if $deg == 0 then null
elif $deg == 1 then -$coefs[0]/$coefs[1]
elif $deg == 2 then first(quadraticRoots)
elif evalPoly(0) == 0 then 0
else 0.001 as $eps
| diffPoly as $deriv
| { x0: $guess, iter: 1, return: null }
| until(.return;
.x0 as $x0
| ($deriv|evalPoly($x0)) as $den
| if $den == 0
then .x0 |= (if . >= 0 then . + $eps else . - $eps end)
else ($coefs|evalPoly($x0)) as $num
| (.x0 - ($num/$den) * $mult) as $x1
| if (($x1 - .x0)|length) <= $tol # abs
then ($x1 | round) as $r
| if (($r - $x1)|length) <= $eps and ($coefs|evalPoly($r)) == 0
then .return = $r
else .return = $x1
end
else .x0 = $x1
end
| if .iter == $maxIter then .return = true
else .iter += 1
end
end)
| if .return != true then .return
else (.x0|round) as $x0
| if ($coefs | evalPoly($x0)) == 0 then $x0
else null
end
end
end;
# Convenience versions of rootPoly/4:
def rootPoly($guess):
rootPoly($guess; 1e-15; 100; 1);
def rootPoly:
rootPoly(0.001; 1e-15; 100; 1);
# Emit a stream of real roots
def roots:
if length==3 then quadraticRoots
else rootPoly as $root
| select($root)
| $root,
# divide by (x - $root)
( divrem( [- $root, 1] ) as [$div, $rem]
| $div
| roots )
end ;
# Given $root is an estimated root of the input polynomial,
# estimate the error dx from deriv = (dY / $root)
# except that if $root or deiv is 0, then use a very small value:
# since 1E-324 == 1E-325 but 1E-323 != 1E-324, we choose 1E-323
def estimatedDeltaX($root):
1E-323 as $tiny
| evalPoly($root) as $dY
| (diffPoly | evalPoly($root)) as $deriv
| if $deriv == 0 then $tiny
else (($dY / $deriv) | length) as $dx
| if $dx == 0 then $tiny
else $dx
end
end ;
def roots_with_errors:
rootPoly as $root
| select($root)
| estimatedDeltaX($root) as $dx
| [$root, $dx],
# divide by (x - $root)
( divrem( [- $root, 1] ) as [$div, $rem]
| $div
| roots_with_errors );
def illustration($text; $poly):
$text,
($poly | roots_with_errors),
"";
"Each root and corresponding estimated error is presented in an array.",
illustration("X^3 = 1"; [-1, 0, 0, 1] ),
illustration("X^3 - 4X^2 - 27X + 90 = 0"; # [-5, 3 ,6]
[90, -27, -4, 1] ),
illustration("(X - 1)^3"; [-1,3,-3,1] ),
illustration("X^2 = 2"; [-2,0,1] ),
illustration("X^3 = 2"; [-2,0,0,1] )
</syntaxhighlight>
{{out}}
<pre>
Each root and corresponding estimated error is presented in an array.
X^3 = 1
[1,1e-323]
X^3 - 4X^2 - 27X + 90 = 0
[3,1e-323]
[-5,1e-323]
[6,1e-323]
(X - 1) ^3
[1,1e-323]
[1,1e-323]
[1,1e-323]
X^2 = 2
[-1.4142135623730951,1.570092458683775e-16]
[1.4142135623730951,1e-323]
X^3 = 2
[1.2599210498948732,1e-323]
</pre>
=={{header|Julia}}==
Julia has extensive matrix math support via the LAPACK library. The Horner function is from the Raku example.
<syntaxhighlight lang = "julia">using LinearAlgebra
""" Cubic polynomial using Horner's method. """
horner(xarr, a, b, c, d) = @. ((a * xarr + b) * xarr + c) * xarr + d
""" Get parameters of characteristic cubic polynomial (with x^3 term 1) from 3x3 matrix """
polynomial3(t) = (1, -tr(t), (tr(t)^2 - tr(t^2)) / 2, -det(t)) # a, b, c, d
# test values
const sin45 = 1 / sqrt(2)
const testmats = [
[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 sin45 -sin45; 0 sin45 sin45],
]
for mat in testmats
display(mat)
a, b, c, d = polynomial3(mat)
print("The characteristic polynomial is: ")
println("x³ ", b < 0 ? "- $(-b)" : " + $b", "x^² ", c < 0 ? "- $(-c)" : "+ $c", "x ",
d < 0 ? "- $(-d)" : "+ $d")
evals = eigvals(mat)
println("The LAPACK library computed eigenvalues are: $evals")
println("Errors are: ", horner(evals, a, b, c, d), "\n")
end
</syntaxhighlight>{{out}}
<pre>
3×3 Matrix{Float64}:
1.0 -1.0 0.0
0.0 1.0 1.0
0.0 0.0 1.0
The characteristic polynomial is: x³ - 3.0x^² + 3.0x - 1.0
The LAPACK library computed eigenvalues are: [1.0, 1.0, 1.0]
Errors are: [0.0, 0.0, 0.0]
3×3 Matrix{Float64}:
-2.0 -4.0 2.0
-2.0 1.0 2.0
4.0 2.0 5.0
The characteristic polynomial is: x³ - 4.0x^² - 27.0x + 90.0
The LAPACK library computed eigenvalues are: [-5.000000000000005, 3.0, 6.0]
Errors are: [-4.831690603168681e-13, 0.0, 0.0]
3×3 Matrix{Float64}:
1.0 -1.0 0.0
0.0 1.0 -1.0
0.0 0.0 1.0
The characteristic polynomial is: x³ - 3.0x^² + 3.0x - 1.0
The LAPACK library computed eigenvalues are: [1.0, 1.0, 1.0]
Errors are: [0.0, 0.0, 0.0]
3×3 Matrix{Float64}:
2.0 0.0 0.0
0.0 -1.0 0.0
0.0 0.0 -1.0
The characteristic polynomial is: x³ + -0.0x^² - 3.0x - 2.0
The LAPACK library computed eigenvalues are: [-1.0, -1.0, 2.0]
Errors are: [0.0, 0.0, 0.0]
3×3 Matrix{Float64}:
2.0 0.0 0.0
0.0 3.0 4.0
0.0 4.0 9.0
The characteristic polynomial is: x³ - 14.0x^² + 35.0x - 22.0
The LAPACK library computed eigenvalues are: [1.0, 2.0, 11.0]
Errors are: [0.0, 0.0, 0.0]
3×3 Matrix{Float64}:
1.0 0.0 0.0
0.0 0.707107 -0.707107
0.0 0.707107 0.707107
The characteristic polynomial is: x³ - 2.414213562373095x^² + 2.414213562373095x - 0.9999999999999998
The LAPACK library computed eigenvalues are: ComplexF64[0.7071067811865475 - 0.7071067811865475im, 0.7071067811865475 + 0.7071067811865475im, 1.0 + 0.0im]
Errors are: ComplexF64[1.1102230246251565e-16 + 1.1102230246251565e-16im, 1.1102230246251565e-16 - 1.1102230246251565e-16im, 2.220446049250313e-16 + 0.0im]
</pre>
=={{header|Phix}}==
=== {{trans|Wren}} ===
Note the six 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. The previously "suspect use of reverse() in eigenvalues(), which kinda brute-forces it to work" is now in Wren too.
<!--(phixonline)-->
<syntaxhighlight lang="phix">
with javascript_semantics
Line 53 ⟶ 396:
function matrix_determinant(sequence m)
// Returns the determinant of the current instance
//
integer l = length(m)
if l=1 then return m[1,1] end if
Line 77 ⟶ 420:
end for
end for
return res
end function
Line 88 ⟶ 425:
function diff_poly(sequence coefs)
// Returns the coefficients of a polynomial following differentiation.
// Polynomials are represented as described in 'eval_poly'
integer c = length(coefs)-1
if c=0 then return {0} end if
Line 139 ⟶ 476:
base -= 1
end while
return res
end function
include complex.e
// assumes real coefficients, highest to lowest order
function quadratic(atom a, b, c)
atom d = b*b - 4*a*c
if d=0 then // single real root
atom root = -b/(a*2)
return {root, root}
elsif d>0 then // two real roots
atom sr = sqrt(d)
d = iff(b<0 ? sr - b : -sr - b)
return {d/(2*a), 2*c/d}
end if
// two complex roots
atom den = 1/(2*a)
complex t1 = complex_new(-b*den, 0),
t2 = complex_new(0, sqrt(-d)*den)
return {complex_add(t1,t2), complex_sub(t1,t2)}
end function
// same assumption as quadratic but 'x' can be real or complex.
function eval_poly(sequence coefs, object x)
object res = coefs[1]
for i=2 to length(coefs) do
if complex(x) then
res = complex_add(complex_mul(res,x),coefs[i])
else
res = res*x + coefs[i]
end if
end for
return res
end function
Line 146 ⟶ 516:
string r = ""
for t=length(si) to 1 by -1 do
if sit!=0 then
if sit=1 and t>1 then
Line 157 ⟶ 527:
sit = abs(sit)
end if
r &= sprintf("%
end if
r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
Line 172 ⟶ 542:
matrix_trace(cofactors(m)),
-matrix_determinant(m)},
roots =
errs =
// find first root
roots
errs
// divide out to get quadratic
sequence den = {-roots[1],1},
// find second and third
errs = append(errs,eval_poly(cp, roots[2]))
string pm = ppf(m,{pp_Nest,1,pp_Indent,8,pp_IntFmt,"%2d"})
return {pm, poly(reverse(cp)), roots, errs}
end function
constant v = 1/sqrt(2)
constant tests = {{{ 1, -1, 0},
{ 0, 1, 1},
Line 200 ⟶ 564:
{{-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}}}
constant fmt = """
Matrix: %s
Line 228 ⟶ 604:
characteristic polynomial: x^3 - 4x^2 - 27x + 90
eigenvalues: {3,6,-5
errors: {0,0,0}
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, 0, 0},
{ 0,-1, 0},
{ 0, 0,-1}}
characteristic polynomial: x^3 - 3x - 2
eigenvalues: {-1,2,-1}
errors: {0,0,0}
Matrix: {{ 2, 0, 0},
{ 0, 3, 4},
{ 0, 4, 9}}
characteristic polynomial: x^3 - 14x^2 + 35x - 22
eigenvalues: {1,11,2}
errors: {0,0,0}
Matrix: {{ 1, 0, 0},
{ 0,0.7071067812,-0.7071067812},
{ 0,0.7071067812,0.7071067812}}
characteristic polynomial: x^3 - 2.41421x^2 + 2.41421x - 1
eigenvalues: {1.0,{0.7071067812,0.7071067812},{0.7071067812,-0.7071067812}}
errors: {2.22044605e-16,{2.22044605e-16,-2.775557561e-16},{2.22044605e-16,2.775557561e-16}}
</pre>
=== {{trans|Raku}} ===
Note "wren2" 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 = {{"wren1", {{ 1, -1, 0}, { 0, 1, 1}, { 0, 0, 1}}},
{"wren2", {{-2, -4, 2}, {-2, 1, 2}, { 4, 2, 5}}},
{"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 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: wren1 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: wren2 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: 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}
4: 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}
5: 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}
6: 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}}==
<small></small>
<syntaxhighlight lang="raku">
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 );
}
</syntaxhighlight>
{{out}}
<pre style="font-size: 12px">
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>]
</pre>
Line 248 ⟶ 921:
{{libheader|Wren-matrix}}
{{libheader|Wren-math}}
{{libheader|Wren-complex}}
{{libheader|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 cubic polynomial after
<syntaxhighlight lang="wren">import "./matrix" for Matrix
import "./math" for Math
import "./complex" for Complex
import "./fmt" for Fmt
Line 292 ⟶ 967:
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)
Line 304 ⟶ 1,009:
cp[3] = -m.det
// find first root
roots
errs
// divide out to get quadratic
var num = Polynom.new(cp[-1..0])
var den = Polynom.new([-roots[0], 1])
// find second
errs.add(evalPoly.call(cp, roots[1]))
return [cp, roots, errs]
}
var
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 (
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:")
Line 366 ⟶ 1,089:
x³ - 4x² - 27x + 90
Its eigenvalues are: [3,
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]
</pre>
| |||