# Newton Raphson- iteration to find roots but COMPLEX as well as REAL..

Ray McAlister piqued my interest in this. I always have to code things myself to really get to grips with them- and LB is ideal for this.

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.

## General situation for finding roots in complex plane.

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.

## Draw the whole complex plane and the fate of starting at every point.

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
'       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    '
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)) ...
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' ...
'   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

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

```