cubic solver
|
09-08-2023, 02:26 PM
(This post was last modified: 03-03-2024 11:58 PM by Albert Chan.)
Post: #1
|
|||
|
|||
cubic solver
This program solve x^3 = a*x + b, for x
Code: cubic_ab(a,b) := Code: cubic(a,b) := Michael Penn's example Let x = sin(θ) + cos(θ), it reduced to cubic: x^3 = 3*x - 11/8 Cas> r := cubic(3, -11/8) Cas> simplify(r) → [1/2, (3*√5-1)/4, (-3*√5-1)/4] Cas> float(Ans) → [0.5, 1.42705098312, −1.92705098312] Cas> float(r) → [1.42705098312, −1.92705098312, 0.5] We have Cas simplify bug here! symbolic and approximate numbers don't match. Luckily, the algorithm doesn't care which cube roots is used (no need for principle root). But, this bug should be fixed. Cas> b := cubic_ab(3, -11/8) [2] → (1/256*(-√15*48*i-176))^(1/3) Cas> polar(float(b)) → 1., −0.77627903074 Cas> polar(simplify(b)) → 1 , atan(√15) Simplify bug had phase angle off by 2*pi/3: Cas> atan(√15) - 2*pi/3.0 → −0.77627903074 Update: cubic_ab(a,b) code explained here |
|||
09-08-2023, 02:44 PM
Post: #2
|
|||
|
|||
RE: cubic solver
Root simplify returned wrong (not principle) root bug existed for XCas 1.5.0
But, it is worse on XCas 1.9.0-31 (mingw32 win32) Result does not match phase *and* magnitude. XCas> b := (1/256*(-√15*48*i-176))^(1/3) XCas> polar(float(b)) 1.0, -0.77627903074 XCas> simplify(polar(b)) // good \(\displaystyle 1\;,\;\frac{-\pi +\arctan \left(\frac{3}{11} \sqrt{15}\right)}{3}\) XCas> simplify(polar(simplify(b))) // bad \(256\;,\;-\pi +\arctan \left(\frac{3}{11} \sqrt{15}\right)\) |
|||
09-08-2023, 06:56 PM
Post: #3
|
|||
|
|||
RE: cubic solver
Yes, there is an additional check in recent versions of Xcas, and here it failed because an argument was not reduced mod 2*pi.
https://github.com/geogebra/giac/commit/...ea8fd94c9d |
|||
06-29-2024, 02:46 PM
(This post was last modified: 07-02-2024 11:06 AM by Albert Chan.)
Post: #4
|
|||
|
|||
RE: cubic solver
Updated cubic solver (x^3=a*x+b), use quadratic solver (x^2=a*x+b) to get alpha^3
Code: quadratic(a,b) := { local d; Code: cubic(a,b) := { (06-29-2024 12:54 AM)Albert Chan Wrote: x^3 - 4*x^2 + 8*x - 8 = 0 // let x = y + 4/3 XCas> cubic(-8/3, 56/27) .+ 4/3. 16/9 [2.0, 1.0+1.73205080757*i, 1.0-1.73205080757*i] Note: (green) cubic discriminant is really quadratic discriminant. Update: both solver now work even if both a, b are 0's |
|||
06-29-2024, 03:29 PM
(This post was last modified: 07-02-2024 11:12 AM by Albert Chan.)
Post: #5
|
|||
|
|||
RE: cubic solver
Previous post was for XCas. Here is for HP Prime Cas
Code: #cas Is there a way to print out discriminant on the same screen? (i.e. without pressing key to get back?) |
|||
07-02-2024, 02:34 PM
Post: #6
|
|||
|
|||
RE: cubic solver
Here is solving cubic with trig identities
cos(3θ) = 4*cos(θ)^3 - 3*cos(θ) sin(3θ) = 3*sin(θ) - 4*sin(θ)^3 Or, to make formula with same pattern 4*cos(θ)^3 = 3*cos(θ) + cos(-3θ) 4*sin(θ)^3 = 3*sin(θ) + sin(-3θ) x^3 = a*x + b // let x = k*y k³*y³ = a*k*y + b // multiply both side by 4/k³, assuming k≠0 4*y³ = 4a/k²*y + 4b/k³ 4*y³ = 3*y + t // where k = √(4a/3), t = 4b/k³ This matched both trig identities. We just need to fix when k=0 (i.e. a=0) For cos version, replace sin/asin with cos/acos. Code: cubic_kt(a,b) := { a:=sqrt(4*a/3); return normal([a, b*4/a^3]) }; x^3 - 4*x^2 + 8*x - 8 = 0 // x = y + 4/3 y^3 = -8/3*y + 56/27 XCas> cubic_sin(-8/3, 56/27) .+ 4/3. [2.0, 1.0+1.73205080757*i, 1.0-1.73205080757*i] |
|||
07-02-2024, 05:12 PM
Post: #7
|
|||
|
|||
RE: cubic solver
(07-02-2024 02:34 PM)Albert Chan Wrote: For cos version, replace sin/asin with cos/acos. For numerical calculations, sin version perhaps is better. (is there exceptions?) But there is a use for cosine, if we want to get real (α, β) from real (k, t) cos(real or purely imaginary) = real y^3 = a*y+b, k = √(4a/3), t = 4b/k³, θ = acos(t)/3 y = k* cos(2n/3*pi - θ) = k*(cos(2n/3*pi)*cos(θ) + sin(2n/3*pi)*sin(θ)) = [cos(2n/3*pi), i*sin(2n/3*pi)] • [k*cos(θ), k*sin(θ)/i] Compare this to radical (α, β) version y = cis(2m/3*pi) * α + cis(-2m/3*pi) * β = [cos(2m/3*pi), i*sin(2m/3*pi)] • [α+β, α−β] Assuming m=n (it doesn't matter if they are not the same), we have: α+β = k*cos(θ) α−β = -i*k*sin(θ) 2α = (α+β) + (α−β) = k*cis(-θ) 2β = (α+β) − (α−β) = k*cis(θ) LHS = one of RHS pairs. (we don't know if m=n, or sort order of (α, β)) k/2*cis(±acos(t)/3) = {α, β} | {α*ω, β/ω} |{α/ω, β*ω} | {β, α} | {β*ω, α/ω} | {β/ω, α*ω} Example, x^3 = x + 1 XCas> a, b := cubic_ab(1.,1.) → [0.337726750973, 0.986991206271] XCas> k, t := cubic_kt(1.,1.) → [1.15470053838, 2.59807621135] XCas> θ := acos(t)/3 → 0.536211995193*i XCas> k/2 * [cis(θ), cis(-θ)] → [0.337726750973, 0.986991206271] |
|||
07-03-2024, 10:17 PM
(This post was last modified: 07-05-2024 09:38 PM by Albert Chan.)
Post: #8
|
|||
|
|||
RE: cubic solver
Previous posts assumed depressed cubic. Here is how to depress a cubic.
This version, if cubic coefficients are integers, so does depressed cubic. Let's define {a,b,c,d} = roots of polynomial, a*x^3 + b*x^2 + c*x + d To force b to 0, we need substitution, x = y - b/A, where A=3a But we don't want fractions. Let's scale the problem. Bonus: cubic coefficient = 1 {a, b, c, d} = {A, 3b, 3c, 3d} = {1, 3b, 3Ac, 3A²d} / A Synthetic Division with offset b will clear 3b coefficient Code: -b> 1 3b 3Ac 3A²d {a,b,c,d} = ({1, 0, C, D} - b) / A Example, super-golden ratio, ψ³ = ψ²+1, we have a, b, c, d = 1, -1, 0, -1 A = 3 * a = 3 C = 3 * (Ac-b²) = 3 * -1 D = 3A²d - b*(C+b²) = -29 α³, β³ = {1, D, -(C/3)³} = {1,-29,1} = (29 ± √(93²-4))/2 = (29 ± 3√93)/2 ψ = ({1, 0, -3, -29} - (-1)) / 3 = (α + β + 1) / 3 Cas> (cubic(3,29) .+ 1) / 3. [1.46557123188, −0.232785615938+0.792551992515*i, −0.232785615938-0.792551992515*i] This matched Michael Penn's ψ solution (08-07-2022 07:55 AM)Thomas Klemm Wrote: Recently I stumbled upon this video by Micheal Penn: What is the super-golden ratio?? Code: function cubic(a,b,c,d) -- a*x^3 + b*x^2 + c*x + d = 0 |
|||
07-03-2024, 11:48 PM
Post: #9
|
|||
|
|||
RE: cubic solver
Based from previous post lua code (depressed cubic code not needed anymore)
HP Prime Cas Code: #cas XCas Code: quadratic(a,b,c) := { local d; |
|||
07-04-2024, 08:37 PM
(This post was last modified: 07-06-2024 03:20 PM by Albert Chan.)
Post: #10
|
|||
|
|||
Quarter Solver
We can use same trick (post #8) to depress quartic, a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
{a, b, c, d, e} = {(A = 4a), 4b, 4c, 4d, 4e} = {1, 4b, 4Ac, 4A²d, 4A³e} / A Synthetic Division to force coefficient 4b to 0 Code: -b> 1 4b 4Ac 4A²d 4A³e {a, b, c, d, e} = ({1, 0, C, D, E} - b) / A y^4 + C*y^2 + D*y + E = 0 (y^2+s)^2 = (2s-C)*y^2 - D*y + (s^2-E) // Add an offset s, making LHS perfect square Let k² = 2s-C ≠ 0, and force RHS as perfect square too. {(k² = 2s-C), -D, (s²-E)} = {1, -D, (D²/4 = (s²-E)*k²)} / k² D² = 4*(s²-E)*k² = ((k²+C)² - 4E) * k² (k²)^3 + 2*C*(k²)^2 + (C²-4E)*(k²) - D² = 0 The reason we solve for k² instead of s is to easily check k²≠0 k² roots all zero ⇒ {C,D,E} all zeroes ⇒ y roots all zero too. (y²+s)² = k² * (y-(D/2)/k²)² // let t = D/2/k (y²+s) = ±(k*y - t) (y^2 - ky + s+t) * (y^2 + ky + s-t) = 0 (s±t), one of them may cause huge cancellation, we only use the good one. The other can be derived by product of all roots = (s+t) * (s-t) = E Last step is simply map y back to x: (x roots) = ((y roots) - b) / A Code: function quartic(a,b,c,d,e) -- a*x^4 + b*x^3 + c*x^2 + d*x + e = 0 lua> quartic(1,2,3,4,5) (0.28781547955764775+1.4160930801719083*I) (-1.287815479557648+0.8578967583284901*I) (0.2878154795576481-1.4160930801719085*I) (-1.287815479557648-0.8578967583284899*I) |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)