but COMPLEX as well as REAL..
If we consider polynomials in x eg f( x) =x^2 -3x +2, we see that it has 'roots' at x =1 and x =2, where f( x) =0.
To find an EXACT root when it happens not to be a nice integer we can draw a nice curve by hand, but much better is to use Newton's method, which basically says make any guess, then a better one will be found by following the tangent at the first guess back to the f( x) =0 axis. This is met in very early calculus class. Basically
x --> x =f( x) /f'( x)
ie replace the guess by its value minus the result of dividing f( x) value at the guess by f'( x) there, where the prime sign ' means first derivative or gradient.
I first coded this a decade or so back in LB.
I've also coded things like Mandelbrots, but Ray brough me back to the Newton Method- but with it generalised to polynomials which have complex roots as well or in place of real ones.
The diagram shows that you can choose ANY point in the complex plane, apply Newton's method, and in the great majority of cases you will converge to one of the possible roots ( occasionally yoou might get an oscillation instead).
The following show the path followed in several cases...
Start near a root- gets there in three steps..
Further from a root- takes seven steps
Zig-zags pretty wildly but gets there...
Another fairly weird path..
Each step is basically a rotation and displacement, applied repeatedly. The further you are from a root the longer it will take to get to the root. And there must be a borderline where a tiny difference in start position could send you to one ot to another root.. hence we bring in fractal behaviour.
If we systematically try all points in the range of interest ( or Monte Carlo fashion, try them all in random order) we can colour code them for which root they ended at. Even better, we can modulate the colour saturation to represent the number of iterations it took.
Unfortunately LB does not have complex numbers and arithmetic operations on them. I wrote some years ago however. And it will be VERY slow. To explore to 2 decimal places all possible starting points between -5 and +5 on real and imaginary axes means 500 *500 starting points, typically iterating 3 to 15 times at each, and with the overhead of my routines, Definitely 'go have a coffee' time, but worth it for the nice graphics. Next improvement is to add the ability to zoom into the fractal areas- they have infinite level of detail!
This page is provisional- code is being developed further and will be added/altered soon. Or PM me...
' **********************************************************************************************************
' ************* Format representing a complex number is a string. "2i-2" is ( 2 -2i) ***************
' ************* ***************
' ************* Supply 'roots' and the polynomial f( x) and derivative f'( x) ***************
' ************* Draw the fate of iterations from all starting points in the mapped area ***************
' ************* ***************
' ************* testRootsInvd7.bas July 2109 tenochtitlanuk ***************
' ************* ***************
' **********************************************************************************************************
' Version to save bursts of typical root seeking iterations...
' to-do's
' add on-screen display of f(x) and routes done
' add save or print screen- save added
' routine to zoom
' routine to plot a dot at start for the given roots done
' move the division to ( corrected) proc 'cDivision('
' trap and debug the divide-by-zero error OK
' 'exit for' if within 'eps' of a root done
' use number of iterations to fade the root colour done
' test with eg f( x) =x^2 -i *x -2 *i, roots 2 and -i, f'( x) =2 *x -i current
' add button to screengrab on demand EDIT to grab an iteration sequence
global p: p =0 ' <<<<<<<<<<<<<<<<<<< for use in drawing iteration route...
allCol$ ="red,green,blue,yellow,lightgray,black,white,lightgray,brown,pink,darkcyan,darkgreen,darkred"
nomainwin
'on error goto [errorHandler]
global delta, type$, m, dist
type$ ="gt-" ' ie graphics, text, both, and with w.g2 suppressed.
firstTime =0
UpperLeftX = 4: UpperLeftY = 4
WindowWidth =1180: WindowHeight =660
textbox #w.tb, 20,580, 500, 30
textbox #w.tb2, 680,580, 500, 30
graphicbox #w.g, 10, 20, 550, 550
graphicbox #w.g2, 602, 20, 550, 550
button #w.b1, "ScrGrab", scr, LR, 590, 20, 80, 30
open "Complex plane- Newton Fractals" for window as #w
#w "trapclose quit"
#w.g "down ; fill 100 120 100"
#w.g2 "down ; fill 100 120 100 ; flush ; color white"
#w "font Courier_New 7"
gosub [grid]
#w.g "flush"
' word$( 1) =f(x), word$( 2) =f'( x), word$( 3) solutions sep'd by '|'.
' f(x), f'(x), root( 1)| root( 2)|...etc
command$ ="(x+2i2)(x-2i-2) OR x^2 -8i, 2x, 2i2|-2i-2"
w$ =word$( command$, 3, ",")
a1$ =word$( w$, 1, "|")
a2$ =word$( w$, 2, "|")
'a3$ =word$( w$, 3, "|")
a1R$ =word$( a1$, 1, "i"): a1I$ =word$( a1$, 2, "i")
a2R$ =word$( a2$, 1, "i"): a2I$ =word$( a2$, 2, "i")
'a3R$ =word$( a3$, 1, "i"): a3I$ =word$( a3$, 2, "i")
R( 1) =val( a1R$): I( 1) =val( a1I$) ' R( n) and I( n) hold real and imaginary parts of roots
R( 2) =val( a2R$): I( 2) =val( a2I$) ' add further roots if polynomial is above quadratic
'R( 3) =val( a3R$): I( 3) =val( a3I$)
delta =0.2
delta =max( 0.02, delta)
eps =0.0001
'while 1
'x =-5 +delta *int( 10 /delta *rnd(1)) ' use for Monte Carlo runs
'y =-5 +delta *int( 10 /delta *rnd(1))
'scan
for x =( 0 -5 -delta /2) to ( 5 +delta /2) step delta
for y =-5 to ( 5 +delta /2) step delta
x$ =str$( x) +"i" +str$( y) ' start at random place in complex ( Real/Imaginary) plane
h$ =x$ ' hold initial start location while x$ iterates..
sR =x: sI =y ' hold previous position on iteration path
gosub [grid2]
'#w.g2 "getbmp scr 1 1 550 550": bmpsave "scr", str$( x) +"_" +str$( y) +"_itern_" +str$(m ) +"_" +str$( time$( "seconds")) +".bmp":wait
'if instr( type$, "-") =0 then
#w.g2 "up ; goto "; 280 +int( sR *50); " "; 280 -int( sI *50)
'end if
for m =0 to 12 ' twelve iterations normally plenty '
topLine$ =cAdd$( cSqu$( x$), "0i-8")
bottomLine$ =cProd$( x$, "2i0")
' multiply top and bottom by complex-conjugate of bottom line to make bottom line real..
CC$ =cConj$( bottomLine$)
topLine$ =cProd$( topLine$, CC$)
bottomLine$ =cProd$( bottomLine$, CC$)
divBy =val( word$( bottomLine$, 1, "i"))
if divBy =0 then exit for
presentR =val( word$( topLine$, 1, "i"))
presentI =val( word$( topLine$, 2, "i"))
R =presentR /divBy
I =presentI /divBy
if abs( R) >1000 then exit for
if abs( I) >1000 then exit for
x$ =cSub$( x$, str$( R) +"i" +str$( I)) ' calculates new approximation
R =val( word$( x$, 1, "i")) ' x --> x -f( x) /f'( x)
I =val( word$( x$, 2, "i"))
if ( abs( R) 0 then #w.g2 "getbmp scr 1 1 550 550": bmpsave "scr", str$( x) +"_" +str$( y) +"_itern_" +str$(m ) +"_" +str$( time$( "seconds")) +".bmp"
next m
'if instr( type$, "t") then #w.tb2 "..hopped to...( "; using( "##.####", sR); ", i";_
' using( "##.####", sI); ")"
'gosub [grid2]
'if p <>0 then #w.g2 "getbmp scr 1 1 550 550": bmpsave "scr", str$( x) +"_" +str$( y) +"_itern" +str$( time$( "seconds")) +".bmp"
if p<>0 then p =0
scan
R =val( word$( x$, 1, "i"))
I =val( word$( x$, 2, "i"))
' find distance of current end point ( R, I) from each root ( R( n), I( n)) ...
' NB add to if >quadratic..<<
dist1 =( ( R -R( 1))^2 +( I -I( 1))^2)^0.5
dist2 =( ( R -R( 2))^2 +( I -I( 2))^2)^0.5
'dist3 =( ( R -R( 3))^2 +( I -I( 3))^2)^0.5
#w.g "backcolor cyan ; color cyan ; down"
m =min( m, 11)
m$ =right$( " " +str$( int( ( 11 -m) /11 *255)), 3)
red$ ="255 " +m$ +" " +m$
green$ =m$ +" 255 " +m$
blue$ =m$ +" " +m$ +" 255"
if dist1 0 and firstTime <> 1 then
#w.g2 "cls ; fill 100 120 100"
#w.g2 "color white" ' subroutine to draw grid and roots
for x2 =-5 to 5
if x2 =0 then #w.g2 "size 4" else #w.g2 "size 1"
#w.g2 "up ; goto "; 280 +x2 *50; " 0 ; down ; goto "; 280 +x2 *50; " 600"
next x2
for y2 =-5 to 5
if y2 =0 then #w.g2 "size 4" else #w.g2 "size 1"
#w.g2 "up ; goto "; 570; " "; 280 +y2 *50; " ; down ; goto 0 "; 280 +y2 *50
next y2
#w.g2 "size 20 ; color yellow ; down"
#w.g2 "set "; 280 -50 *2; " "; 280 +50 *2
#w.g2 "set "; 280 +50 *2; " "; 280 -50 *2
'#w.g2 "set "; 280 +50 *0; " "; 280 +50 *1
#w.g2 "size 1"
'end if
firstTime =1
return
[grid]
#w.g "color white" ' subroutine to draw grid and roots
for x =-5 to 5
if x =0 then #w.g "size 4" else #w.g "size 1"
#w.g "up ; goto "; 280 +x *50; " 0 ; down ; goto "; 280 +x *50; " 600"
next x
for y =-5 to 5
if y =0 then #w.g "size 4" else #w.g "size 1"
#w.g "up ; goto "; 570; " "; 280 +y *50; " ; down ; goto 0 "; 280 +y *50
next y
#w.g "size 20 ; color yellow ; down"
if instr( type$, "-") =0 then
#w.g2 "set "; 280 -50 *2; " "; 280 +50 *2
#w.g2 "set "; 280 +50 *2; " "; 280 -50 *2
'#w.g2 "set "; 280 +50 *0; " "; 280 +50 *1
end if
#w.g "size 1"
return
sub quit h$
close #w
end
end sub
end
sub scr h$
'#w.g "up ; goto 120 20 ; down ; color green ; backcolor 100 120 100 ; font 7"
'#w.g "\"; command$
'#w.g "flush"
'#w.g "getbmp scr 1 1 1200 650"
'bmpsave "scr", "test4.bmp"
p =1
end sub
[errorHandler]
print "Error string is " + chr$(34) + Err$ + chr$(34)
print "Error number is ";Err
end
' Collection of functions of a complex number a$ expressed in 'real i imag' ...
' cAdd$( a$, b$)
' cSub$( a$, b$)
' cConj$( a$)
' cProd$( a$, b$)
' cSqu$( a$)
' cCube$( a$)
' check( a$)
' cDiv( a$, b$)
' cPrint(
' a$ uses global selector "t"|"g"|"gt" for text only, graphics only or both
function i$( k$)
R =val( word$( k$, 1, "i"))
I =val( word$( k$, 2, "i"))
tmp =I
I =R
R =0 -tmp
i$ =str$( R) +"i" +str$( I)
end function
function cConj$( I$)
R =val( word$( I$, 1, "i"))
I =val( word$( I$, 2, "i"))
I =0 -I
cConj$ =str$( R) +"i" +str$( I)
end function
function cSqu$( i$)
cSqu$ =cProd$( i$, i$)
end function
function cCube$( i$)
cCube$ =cProd$( i$, cSqu$( i$))
end function
function check( c$)
check =0
sepPosn =instr( c$, "i")
select case
case sepPosn =0
print " Malformed c number! No 'i'.": check =1
case str$( val( word$( c$, 1, "i"))) <>word$( c$, 1, "i")
print " Malformed c number! Number error in real part.": check =2
case str$( val( word$( c$, 2, "i"))) <>word$( c$, 2, "i")
print " Malformed c number! Number error in imaginary part.": check =3
end select
end function
function cDiv$( a$, b$) ' <<<<< TBC <<<<<<<<<<<<<<<<
' multiply top and bottom by c conjugate of denominator ( denominator becomes a pure real) and calculate.
aR =val( word$( a$, 1, "i")): aI =val( word$( a$, 2, "i"))
CC$ =cConj$( b$)
CCR =val( word$( CC$, 1, "i")): CCI =val( word$( CC$, 2, "i"))
if ( abs( CCR) <1e6) and ( abs( CCI) >1e-6) then
R =aR /CCR
I =aI /CCR
else
R =1E6
Y =1E6
end if
cDiv$ =str$( R) +"i" +str$( I)
end function
function cProd$( a$, b$)
aR =val( word$( a$, 1, "i")): aI =val( word$( a$, 2, "i"))
bR =val( word$( b$, 1, "i")): bI =val( word$( b$, 2, "i"))
R =aR *bR -aI *bI
I =aR *bI +aI *bR
cProd$ =str$( R) +"i" +str$( I)
end function
function cAdd$( a$, b$)
aR =val( word$( a$, 1, "i")): aI =val( word$( a$, 2, "i"))
bR =val( word$( b$, 1, "i")): bI =val( word$( b$, 2, "i"))
R =aR +bR
I =aI +bI
cAdd$ =str$( R) +"i" +str$( I)
end function
function cSub$( a$, b$)
aR =val( word$( a$, 1, "i")): aI =val( word$( a$, 2, "i"))
bR =val( word$( b$, 1, "i")): bI =val( word$( b$, 2, "i"))
R =aR -bR
I =aI -bI
cSub$ =str$( R) +"i" +str$( I)
end function