Roots of a cubic polynomial: Difference between revisions
Content added Content deleted
imported>Raygis (add example) |
(Corrected first term of matrix example and added Wren solution.) |
||
Line 8: | Line 8: | ||
a*x**3 + b*x**2 + c*x + d = 0 |
a*x**3 + b*x**2 + c*x + d = 0 |
||
matrix: [[1,- |
matrix: [[1,-1,0], [0,1,-1],[0,0,1]] |
||
polynomial: [a,b,c,d] =[1, -3,+3, -1] |
polynomial: [a,b,c,d] =[1, -3,+3, -1] |
||
roots: [1,1,1] |
roots: [1,1,1] |
||
=={{header|Wren}}== |
|||
{{libheader|Wren-matrix}} |
|||
{{libheader|Wren-math}} |
|||
{{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 polynomial after each root is found. |
|||
<syntaxhighlight lang="wren">import "./matrix" for Matrix |
|||
import "./math" for Math |
|||
import "./fmt" for Fmt |
|||
class Polynom { |
|||
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]) |
|||
} |
|||
toString { "Polynomial(%(_factors.join(", ")))" } |
|||
} |
|||
var eigenvalues = Fn.new { |m| |
|||
var roots = List.filled(3, 0) |
|||
var errs = List.filled(3, 0) |
|||
// find the characteristic polynomial |
|||
var cof = m.cofactors |
|||
var cp = List.filled(4, 0) |
|||
cp[0] = 1 |
|||
cp[1] = -m.trace |
|||
cp[2] = cof[0, 0] + cof[1, 1] + cof[2, 2] |
|||
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 [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("\nIts eigenvalues are: $n", eigs[0]) |
|||
Fmt.print("and the corresponding errors are: $n\n", eigs[1]) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
For matrix: |
|||
| 1 -1 0| |
|||
| 0 1 1| |
|||
| 0 0 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| |
|||
Its eigenvalues are: [3, -5, 6] |
|||
and the corresponding errors are: [0, 0, 0] |
|||
</pre> |
Revision as of 19:39, 23 November 2023
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]
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 {
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])
}
toString { "Polynomial(%(_factors.join(", ")))" }
}
var eigenvalues = Fn.new { |m|
var roots = List.filled(3, 0)
var errs = List.filled(3, 0)
// find the characteristic polynomial
var cof = m.cofactors
var cp = List.filled(4, 0)
cp[0] = 1
cp[1] = -m.trace
cp[2] = cof[0, 0] + cof[1, 1] + cof[2, 2]
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 [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("\nIts eigenvalues are: $n", eigs[0])
Fmt.print("and the corresponding errors are: $n\n", eigs[1])
}
- Output:
For matrix: | 1 -1 0| | 0 1 1| | 0 0 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| Its eigenvalues are: [3, -5, 6] and the corresponding errors are: [0, 0, 0]