Attempt at solving roots of polynomial

Recommended Posts

Here is my attempt at writing a script that will solve any roots. It uses the Aberth method

http://en.wikipedia.org/wiki/Aberth_method And it can even solve complex roots. It only goes up to the third degree for right now but I will improve it.

Here is the code

```#include <GUIConstantsEx.au3>
#include "complex.au3"

\$gui = GUIcreate("Root Solver", 900, 500)
\$exit = GUICtrlCreateButton("EXIT", 100, 100, 50, 50)
\$input = GUICtrlCreateinput("",20, 50, 220, 30)
\$ok = GUICtrlCreatebutton("OK", 30, 100, 50, 50)
\$list = GUICtrlCreateList("", 280, 50, 380, 400)
\$list2 = GUICtrlCreateList("",680, 50, 180, 400)
GUICtrlCreateLabel("Iteration Trail",680, 25)
GUICtrlCreateLabel("Put in standard form, put 0x^2 if there is no x^2",15, 22)

GUIsetstate(@SW_SHOW)

While 1
\$msg = GUIgetmsg()
Select
Case \$msg = \$exit
Exit
Case \$msg = \$GUI_EVENT_CLOSE
Exit
Case \$msg = \$ok
GUICtrlSetData(\$list, "")
GUICtrlSetData(\$list2, "")
polynomial()
Endselect
WEnd

Func polynomial()
\$array = Stringsplit(\$userinput, "x")
Select
Case StringinStr(\$array[2], "^2") = 1
second()
Case StringinStr(\$array[2], "^3") = 1
third()
EndSelect
EndFunc

Func second()
\$array = Stringsplit(\$userinput, "x")
For \$i = 1 to \$array[0]
\$array[\$i] = Stringstripws(\$array[\$i], 8)
next

If Stringinstr(\$array[1], "-") = 1 Then
\$array[1] = (-1)*(StringTrimleft(\$array[1], 1))
If StringTrimLeft(\$array[1], 1) = "0" Then
\$array[1] = -1
EndIf
EndIf

If StringLen(\$array[1]) = 0 Then
\$array[1] = 1
EndIf

If StringInStr(\$array[2], "-") > 0 Then
\$b = (-1)*(StringTrimLeft(\$array[2], 3))
If StringTrimleft(\$b, 1) = "0" Then
\$b = -1
EndIf
ElseIf StringinStr(\$array[2], "+") > 0 Then
\$b = StringTrimLeft(\$array[2], 3)
EndIf

If Stringlen(\$b) = 0 Then
\$b = 1
EndIf

If StringinStr(\$array[3], "+") = 1 Then
\$c = StringtrimLeft(\$array[3], 1) + 0.0
ElseIf StringinStr(\$array[3], "-") = 1 Then
\$c = (-1)*(StringTrimLeft(\$array[3], 1))
EndIf

If StringTrimLeft(\$array[3], 1) = 0 Then
\$c = \$array[3]
EndIf

\$a = complex(\$array[1], 0)
\$b = complex(\$b, 0)
\$c = complex(\$c, 0)
\$two = complex(2,0)
\$l = complex(4, 0)
\$i = complex(-1,0)
\$g = complex_pow(\$b, 2)
\$l = complex_mul(complex_mul(\$a, \$c), \$l)
\$g = complex_pow(complex_minus(\$g, \$l), 0.5)
\$b = complex_mul(\$b, \$i)
\$x1 = complex_div(complex_add(\$b, \$g), complex_mul(\$a, \$two))
\$x2 = complex_div(complex_minus(\$b, \$g), complex_mul(\$a, \$two))
GUICtrlSetData(\$list, "x = "&\$x1)
GUICtrlSetData(\$list, "x = "&\$x2)
EndFunc

;uses Aberth Method, it is much better than Newton or Durand-Kerner
Func third()
\$array = Stringsplit(\$userinput, "x")
For \$i = 1 to \$array[0]
\$array[\$i] = Stringstripws(\$array[\$i], 8)
next

For \$i = 2 to \$array[0] - 1
If StringInStr(\$array[\$i], "-") > 0 Then
\$array[\$i] = (-1)*(StringTrimLeft(\$array[\$i], 3))
If StringTrimleft(\$array[\$i], 1) = "0" Then
\$array[\$i] = -1
EndIf
ElseIf StringinStr(\$array[\$i], "+") > 0 Then
\$array[\$i] = StringTrimLeft(\$array[\$i], 3)
EndIf
Next

For \$i = 2 to \$array[0] - 1
If Stringlen(\$array[\$i]) = 0 Then
\$array[\$i] = 1
EndIf
Next

For \$i = 1 to \$array[0]
If StringLen(\$array[\$i]) = 0 Then
\$array[\$i] = 1
EndIf
Next

For \$i = 1 to \$array[0] step 3
If StringInStr(\$array[\$i], "-") > 0 Then
\$array[\$i] = (-1)*(StringTrimLeft(\$array[\$i], 1))
If StringTrimLeft(\$array[\$i], 1) = "0" Then
\$array[\$i] = -1
EndIf
EndIf
Next

\$z = complex(Random(), Random()); Must get better seeds
\$z1 = complex(Random(), Random())
\$z2 = complex(Random(), Random())
;Cannot change the coefficients to negative
\$a = complex(\$array[1], 0)
\$b = complex(\$array[2], 0)
\$c = complex(\$array[3], 0)
\$d = complex(\$array[4], 0)
\$one = complex(1, 0)
\$l = complex(-1, 0)
\$i = 1
While abs(complex_extract_re(complex_add(complex_mul(\$a, complex_pow(\$z, 3)),complex_mul(\$b, complex_pow(\$z, 2)),complex_mul(\$c, \$z), \$d))) > 1e-10
\$v = complex_add(complex_div(\$one, complex_minus(\$z, \$z1)), complex_div(\$one, complex_minus(\$z, \$z2)))
\$m = complex_div(complex_add(complex_mul(\$a, complex_pow(\$z, 3)),complex_mul(\$b, complex_pow(\$z, 2)),complex_mul(\$c, \$z), \$d) , complex_add(complex_mul(complex(3, 0),\$a, complex_pow(\$z, 2)), complex_mul(complex(2, 0),\$b,\$z), \$c))
\$u = complex_mul(\$l, complex_div(\$m, complex_minus(\$one, complex_mul(\$m, \$v))))
\$z = complex_add(\$z, \$u);
\$v = complex_add(complex_div(\$one, complex_minus(\$z1, \$z)), complex_div(\$one, complex_minus(\$z1, \$z2)))
\$m = complex_div(complex_add(complex_mul(\$a, complex_pow(\$z1, 3)),complex_mul(\$b, complex_pow(\$z1, 2)),complex_mul(\$c, \$z1), \$d) , complex_add(complex_mul(complex(3, 0),\$a, complex_pow(\$z1, 2)), complex_mul(complex(2, 0),\$b,\$z1), \$c))
\$u = complex_mul(\$l, complex_div(\$m, complex_minus(\$one, complex_mul(\$m, \$v))))
\$z1 = complex_add(\$z1, \$u);
\$v = complex_add(complex_div(\$one, complex_minus(\$z2, \$z1)), complex_div(\$one, complex_minus(\$z2, \$z)))
\$m = complex_div(complex_add(complex_mul(\$a, complex_pow(\$z2, 3)),complex_mul(\$b, complex_pow(\$z2, 2)),complex_mul(\$c, \$z2), \$d) , complex_add(complex_mul(complex(3, 0),\$a, complex_pow(\$z2, 2)), complex_mul(complex(2, 0),\$b,\$z2), \$c))
\$u = complex_mul(\$l, complex_div(\$m, complex_minus(\$one, complex_mul(\$m, \$v))))
\$z2 = complex_add(\$z2, \$u)
GUICtrlSetData(\$list2, "Loop "&\$i&" = "&\$z&", "&\$z1&", "&\$z2)
\$i = \$i + 1
If abs(complex_extract_re(complex_add(complex_mul(\$a, complex_pow(\$z, 3)),complex_mul(\$b, complex_pow(\$z, 2)),complex_mul(\$c, \$z), \$d))) < 1e-10 or abs(complex_extract_re(complex_add(complex_mul(\$a, complex_pow(\$z1, 3)),complex_mul(\$b, complex_pow(\$z1, 2)),complex_mul(\$c, \$z1), \$d))) < 1e-10 or abs(complex_extract_re(complex_add(complex_mul(\$a, complex_pow(\$z2, 3)),complex_mul(\$b, complex_pow(\$z2, 2)),complex_mul(\$c, \$z2), \$d))) < 1e-10 Then
ExitLoop
EndIf
Wend

GUICtrlSetData(\$list, "x1 = "&\$z)
GUICtrlSetData(\$list, "x2 = "&\$z1)
GUICtrlSetData(\$list, "x3 = "&\$z2)
EndFunc```

Here is the code for the complex number include file. It is my own UDF that I created because I couldnt find any complex numbers UDFs.

```#include <Math.au3>
\$pi = 3.141592653589793

Global const \$ni = complex(0,0)
Global const \$gi = complex(1,0)

Func complex(\$re, \$im)
\$complex = ""&\$re&","&\$im&"i"
return \$complex
EndFunc

;returns the magnitude and the argument(theta). This will be used a lot in conjunction with other
;functions
Func complex_polar(\$x)
\$x = StringSplit(\$x, ","); Splits the complex to a, bi
\$x[2] = StringTrimright(\$x[2], 1); Strips the i off of the imaginary number
\$x[1] = \$x[1] + 0.0; needs to add 0.0 to make it a float for squaring, type conversion
\$x[2] = \$x[2] + 0.0
\$mag = sqrt(\$x[1]^2 + \$x[2]^2)
\$theta = (atan(abs(\$x[2])/abs(\$x[1])))*(180/\$pi)
;this is to select the quadrant for the magnitude and theta. 3 -7i would be in the 4th quadrant
Select
Case \$x[1] > 0 and \$x[2] >= 0
\$theta = \$theta
Case \$x[1] < 0 and \$x[2] >= 0
\$theta = 180 - \$theta
Case \$x[1] < 0 and \$x[2] <= 0
\$theta = 180 + \$theta
Case \$x[1] > 0 and \$x[2] <= 0
\$theta = 360 - \$theta
EndSelect
If \$theta <> \$theta Then
\$theta = 0
EndIf
\$complex = ""&\$mag&","&\$theta&""
Return \$complex
EndFunc

;when this function is used with polar, you can do a lot of things
Func complex_rect(\$complex)
\$complex = StringSplit(\$complex, ",")
\$real = Round(\$real, 12); rounds the number, can change if you want more precision
\$imag = Round(\$imag, 12)
If abs(\$real) < 1e-10 Then; more ifs and thens for more precision handling >_>
\$real = 0
EndIf
If abs(\$imag) < 1e-10 Then
\$imag = 0
EndIf
\$complex = ""&\$real&","&\$imag&"i"
Return \$complex
EndFunc

;this is an example of using the function polar and rect together.
;this can raise to a complex number to any power, even to float like 2.3
Func complex_pow(\$complex, \$pow)
\$complex = complex_polar(\$complex)
\$complex = StringSplit(\$complex, ",")
\$complex[1] = \$complex[1]^\$pow
\$complex[2] = \$complex[2]*\$pow
\$complex = ""&\$complex[1]&","&\$complex[2]&""
\$complex = complex_rect(\$complex)
Return \$complex
EndFunc

Func complex_conj(\$complex1)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$complex1[2] = (-1)*\$complex1[2] + 0.0
\$real = \$complex1[1]
\$imag = \$complex1[2]
\$complex = ""&\$real&","&\$imag&"i"
Return \$complex
EndFunc

Func complex_add(\$complex1, \$complex2, \$complex3 = \$ni, \$complex4 = \$ni, \$complex5 = \$ni)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$complex2 = StringSplit(\$complex2, ",")
\$complex2[2] = StringTrimRight(\$complex2[2], 1)
\$complex3 = StringSplit(\$complex3, ",")
\$complex3[2] = StringTrimRight(\$complex3[2], 1)
\$complex4 = StringSplit(\$complex4, ",")
\$complex4[2] = StringTrimRight(\$complex4[2], 1)
\$complex5 = StringSplit(\$complex5, ",")
\$complex5[2] = StringTrimRight(\$complex5[2], 1)
\$real = \$complex1[1] + \$complex2[1] + \$complex3[1] + \$complex4[1] + \$complex5[1]
\$imag = \$complex1[2] + \$complex2[2] + \$complex3[2] + \$complex4[2] + \$complex5[2]
\$complex = ""&\$real&","&\$imag&"i"
Return \$complex
EndFunc

Func complex_minus(\$complex1, \$complex2, \$complex3 = \$ni, \$complex4 = \$ni, \$complex5 = \$ni)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$complex2 = StringSplit(\$complex2, ",")
\$complex2[2] = StringTrimRight(\$complex2[2], 1)
\$complex3 = StringSplit(\$complex3, ",")
\$complex3[2] = StringTrimRight(\$complex3[2], 1)
\$complex4 = StringSplit(\$complex4, ",")
\$complex4[2] = StringTrimRight(\$complex4[2], 1)
\$complex5 = StringSplit(\$complex5, ",")
\$complex5[2] = StringTrimRight(\$complex5[2], 1)
\$real = \$complex1[1] - \$complex2[1] - \$complex3[1] - \$complex4[1] - \$complex5[1]
\$imag = \$complex1[2] - \$complex2[2] - \$complex3[2] - \$complex4[2] - \$complex5[2]
\$complex = ""&\$real&","&\$imag&"i"
Return \$complex
EndFunc

Func complex_mul(\$complex1, \$complex2, \$complex3 = \$gi, \$complex4 = \$gi, \$complex5 = \$gi)
\$complex1 = complex_polar(\$complex1)
\$complex2 = complex_polar(\$complex2)
\$complex3 = complex_polar(\$complex3)
\$complex4 = complex_polar(\$complex4)
\$complex5 = complex_polar(\$complex5)
\$complex1 = StringSplit(\$complex1, ",")
\$complex2 = StringSplit(\$complex2, ",")
\$complex3 = StringSplit(\$complex3, ",")
\$complex4 = StringSplit(\$complex4, ",")
\$complex5 = StringSplit(\$complex5, ",")
\$mag = \$complex1[1]*\$complex2[1]*\$complex3[1]*\$complex4[1]*\$complex5[1]
\$theta = \$complex1[2] + \$complex2[2] + \$complex3[2] + \$complex4[2] + \$complex5[2]
\$complex = complex_rect(""&\$mag&","&\$theta&"")
Return \$complex
EndFunc

Func complex_div(\$complex1, \$complex2)
\$complex3 = complex_conj(\$complex2)
\$complex1 = complex_mul(\$complex1, \$complex3)
\$complex2 = complex_mul(\$complex2, \$complex3)
\$complex1 = StringSplit(\$complex1, ",")
\$complex2 = StringSplit(\$complex2, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$real = \$complex1[1]/\$complex2[1]
\$imag = \$complex1[2]/\$complex2[1]
\$complex = ""&\$real&","&\$imag&"i"
Return \$complex
EndFunc

Func complex_abs(\$complex1)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$complex1[1] = \$complex1[1] + 0.0
\$complex1[2] = \$complex1[2] + 0.0
\$mag = sqrt(\$complex1[1]^2 + \$complex1[2]^2)
Return \$mag
EndFunc

Func complex_extract_re(\$complex1)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[1] = \$complex1[1] + 0.0
Return \$complex1[1]
EndFunc

Func complex_extract_im(\$complex1)
\$complex1 = StringSplit(\$complex1, ",")
\$complex1[2] = StringTrimRight(\$complex1[2], 1)
\$complex1[2] = \$complex1[2] + 0.0
Return \$complex1[2]
EndFunc```

Dont be too harsh, it is my first attempt at writing complex numbers UDF. You need to write x^3 -2 as x^3 +0x^2 +0x -2. Btw can anyone show me the optimizing process or equation for choosing the seeds for the Aberth method? Right now the seeds are chosen randomly. Thanks

Edited by Ultima2

Share on other sites

So as a math major who's taken numerical analysis I enjoyed seeing this-- not much non-trivial math appears on these forums at all. That said, I feel like you'd have an easier time of this if you either (1) Used a language that inherently understood imaginary numbers (i.e. matlab) or (2) Used an object oriented language where you could create Polynomial and Complex classes. It's not that I think that what your doing is at all bad, just that your using AutoIT for something that it really wasn't designed for. You wouldn't use a tack hammer to drive in a 16d nail, so I wouldn't suggest using AutoIT for doing complicated numerical work.

That said I will PM you the original AMS article this approximation was published in. And I note here that this scheme is cubically convergent... which is pretty damn good. Makes me wonder why we didn't cover it in class.

Edited by Wus

Thank you, Wus

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account