Hi everybody

now let us have a taste of complex numbers.

They are powerful tools in mathematics, and they are used in this short program to solve a third degree equation (or cubic equation).

Program is tested, but nobody can be sure 100%

By By

**Code:**'Cubic equation solver
'by prof. Daniele Mazza, mazzad50@gmail.com or www.molecularmodels.eu
'Feel free to contact me
'
'This software is presented as is, with no
'warrantees expressed or implied.
'
'This software is FREEWARE. You may use it
'any way you want, including in commercial
'applications. A credit to the author is
'not necessary, but can be included if you like.
'
'
'
' *** CUBIC EQUATION SOLVER ***
'
'
' Equation is in its general form a*x^3 + b*x^2 + c*x + d = 0
' It is solved by using the famous CARDANO's formula
' _________________________ _________________________
' \ 3 / ___________ \ 3 / ___________
' \ / q \ /q^2 p^3 \ / q \ /q^2 p^3
' x = \ / - --- + \ / --- + ---- + \ / - --- - \ / --- + ---- - b/(3*a)
' \/ 2 \/ 4 27 \/ 2 \/ 4 27
'
' where p = -(b^2/(3*a^2) + c/a q = 2*b^3/(27*a^3) - b*c/(3*a^2) + d/a
' and the discriminant 'D' is equal to 'q^2/4 + p^3/27'
'
' This formula must make use of complex numbers, in order to find all 3 solutions !!
'
'
pi = 3.1415926535
ti = 180/pi
dim n(6,2), m(3,2)' imaginary numbers, their real part is n(',1) imaginary n(',2)
dim th(6) ' corresponding 'phi' angles in complex plane angles (radians)
print " **** CUBIC EQUATION SOLVER ****":print
print " enter coefficients for a*x^3 + b*x^2 + c*x + d = 0"
input " a ";a
input " b ";b
input " c ";c
input " d ";d
p = c/a - b^2/(3*a^2)
q = 2*b^3/(27*a^3) - b*c/(3*a^2) + d/a
disc = q^2/4 + p^3/27
select case
' ----------------------------- D < 0 (3 real solutions)
case disc<0
a1 = -1*q/2 ' real part
b1 = sqr(-1*(disc)) ' imaginary part
m1 = sqr(a1^2 + b1^2)' modulus of complex number a1 + i*b1 where i = sqr(-1)
m2 = m1^(1/3) ' modulus after extraction of cubic root
if a1<0 then add = pi else add = 0 ' correction needed to bring angles in range 0° - 180° (anticlockwise)
if a1 = 0 then ' just in case, to avoid atn(infinite)
th1 = pi/2
else
th1 = atn(b1/a1) + add ' th1 is 'phi' angle in the complex plane of the argument of the first cubic root
' whils the 'phi' angle of the argument of the second cubic root is -th1
end if
' extraction of first cubic root (three values)
for i = 1 to 3
th(i) = th1/3 + (i-1)*2/3*pi ' dividing th1 by 3 one must care of 120° (360°/3) multiples of it
n(i,1) = cos(th(i))*m2 ' real part
n(i,2) = sin(th(i))*m2 ' imaginary part
next i
' extraction of second cubic root (three values)
for i = 4 to 6
th(i) = -1*th1/3 + (i-4)*2/3*pi
n(i,1) = cos(th(i))*m2 ' real part
n(i,2) = sin(th(i))*m2 ' imaginary part
next i
' among 9 (3 x 3) combinations one checks if imaginary parts are complementary and therefore valid solutions
k = 0:for i = 1 to 3
for j = 4 to 6
if abs ((th(i)+th(j)) mod (2*pi)) < 0.00001 then
k = k + 1
m(k,1) = n(i,1) + n(j,1) - b/(3*a)
m(k,2) = n(i,2) + n(j,2)
end if
next j
next i
' ----------------------------- D > 0 (1 real solution and 2 imaginary)
' --------------------------- or D = 0 (3 real solutions of which at least 2 are coincident)
case disc>=0
r1 = -1*q/2 + sqr(disc)
if r1>0 then
q1 = r1^(1/3)
else
q1 = -1*(-1*r1)^(1/3)
end if
r2 = -1*q/2 - sqr(disc)
if r2>0 then
q2 = r2^(1/3)
else
q2 = -1*(-1*r2)^(1/3)
end if
s1 = sqr(3)/2
m(1,1) = q1 + q2 - b/(3*a) : m(1,2) = 0
m(2,1) = -1*(q1 + q2)/2 - b/(3*a) : m(2,2) = s1*(q1 - q2)
m(3,1) = -1*(q1 + q2)/2 - b/(3*a) : m(3,2) = s1*(q2 - q1)
end select
' print solutions
print:print " The three solution are :"
print " +----------------+----------------+"
print " | real part | imaginary part |"
print " +----------------+----------------+"
for i = 1 to 3
print " | ";using("###.#####",m(i,1));tab(19);"| ";using("###.#####",m(i,2));tab(36);"|"
print " +----------------+----------------+"
next i
end