Quote: "If an 8 year old child can understand the answers to the questions; then I think I'd also understand...
"
If that's in reference to my soon-to-be attempt at teaching, it's 8th grader and this definitely wouldn't be in there!
I'm planning on writing tutorials for this stuff... but I've been really lazy in revamping my website.
Basically, newton's method goes like this. At any point on a continuous (and differentiable) curve, you have a point on the line (at (x,f(x)) ), and a line tangent to that point with slope f'(x), where f' is the derivative of f.
The tangent line intersects the x-axis, and this is a good estimation of the root of the function.
Wikipedia has a good animated gif on this.
Given a slope and a point, the equation of a line through the given point is:
y=m(x-x0)+y0
where (x0,y0) is the point.
So where y=0 we have
0=f'(x) *(xnew-x)+f(x)
solving for xnew (the point where y=0 approximates the point where f(x)=0)
-f'(x)*xnew=-x*f'(x)+f(x)
xnew=x-f(x)/f'(x)
So if we iterate this process again and again, we have something like:
x1=x0-f(x0)/f'(x0)
x2=x1-f(x1)/f'(x1)
x3=x1-f(x2)/f'(x2)
etc.
and x100 will be a really good approximation of the root.
The is an actual proof for newton's method converging, but this explanation is best. We take the linear approximation of the line using the slope at a point, and the root of the line is approximately the root of the function.
As you might imagine, there are plenty of places where this method won't converge. If there is no root, it will jump around forever. If there is no root closeby, it might jump around forever. It might cycle (so after n steps its back where it started), or it might do other weird stuff.
Cycling? Nonconvergence? This starts to sound a lot like the mandelbrot set. In fact, if you plug in a complex polynomial and complex starting values, you get a fractal!
The logarithms in the root finding code? I have no idea why that works. It's a trick I found online to get the shading smooth.
Also, the slider bar is not used in this, I just forgot to remove the code xD
All the unitleft/unitright values are just for getting the unit space into screen/pixel space.
Here's a code snippet that allows you to observe the process:
`%Project Title%
`%Source File Name%
`======================
global w as integer
global h as integer
w=screen width()
h=screen height()
type complex
a as float
b as float
endtype
global retval as complex
dim cpoly(0) as complex
dim iterarray(0) as complex
global numiterations as integer
global plotw as integer
global ploth as integer
global plotoffset as integer
global unittop as float
global unitbottom as float
global unitleft as float
global unitright as float
unittop=2.0
unitbottom=-1.0
unitleft=-2.0
unitright=4.0
ploth=h-300
numiterations=100
plotw=((unitright-unitleft)*ploth)/(unittop-unitbottom) //preserve aspect ratio
plotoffset=(w-plotw)/2
sync on
sync rate 0
setRandomPolynomial(4,5)
resetArray()
do
cls
iterateArray()
if spacekey()
setRandomPolynomial(rnd(4)+2,rnd(10)+1)
resetArray()
endif
lock pixels
drawArray()
unlock pixels
wait 500
sync
loop
function drawArray()
sz=array count(iterarray())
for n=0 to sz
x=xUnitToPixel(iterarray(n).a)
y=yUnitToPixel(iterarray(n).b)
dot x,y
next n
endfunction
function iterateArray()
sz=array count(iterarray())
for n=0 to sz
mytmp as complex
mytmp=iterarray(n)
iterateNewtonsMethod(mytmp)
iterarray(n).a=retval.a
iterarray(n).b=retval.b
next n
endfunction
function resetArray()
undim iterarray()
dim iterarray(-1) as complex
num=0
for x=w/2-200 to w/2+200 step 5
for y=h/2-200 to h/2+200 step 5
array insert at bottom iterarray()
iterarray(num).a=xpixeltounit(x)
iterarray(num).b=ypixeltounit(y)
inc num,1
next y
next x
endfunction
function xpixeltounit(x as float)
v#=(x-screen width()/2)*8.0/screen width()
endfunction v#
function ypixeltounit(x as float)
v#=-(x-screen height()/2)*8.0/screen width()
endfunction v#
function xunittopixel(x as float)
v#=x*screen width()/8+screen width()/2
endfunction v#
function yunittopixel(x as float)
v#=-screen width()*x/8+screen height()/2
endfunction v#
function setclassicpolynomial() //x^3+1
undim cpoly()
dim cpoly(3)
cpoly(0).a=1
cpoly(0).b=0
cpoly(1).a=0
cpoly(1).b=0
cpoly(2).a=0
cpoly(2).b=0
cpoly(3).a=1
cpoly(3).b=0
endfunction
function setRandomPolynomial(degree,coefficientsize)
undim cpoly()
dim cpoly(degree) as complex
for n=0 to degree
cpoly(n).b=rnd(coefficientsize*2)-coefficientsize
cpoly(n).a=rnd(coefficientsize*2)-coefficientsize
next n
endfunction
function iterateNewtonsMethod(pt as complex)
f as complex //represents f(pt)
d as complex //represents f'(pt)
p as complex //used to represent pt then pt^2 then pt^3 etc.
//r as complex r used to stand for "new root", but we don't need it since we only calculate 1 iteration
//slightly more bare bones, 1 iteration only
f.a=0
f.b=0
d.a=0
d.b=0
p.a=1
p.b=0
for j=0 to array count(cpoly())-1
//f=f+p*cpoly(j)
f.a=f.a+p.a*cpoly(j).a-p.b*cpoly(j).b
f.b=f.b+p.b*cpoly(j).a+p.a*cpoly(j).b
//d=d+p*cpoly(j+1)*(j+1)
d.a=d.a+(p.a*cpoly(j+1).a-p.b*cpoly(j+1).b)*(j+1)
d.b=d.b+(p.b*cpoly(j+1).a+p.a*cpoly(j+1).b)*(j+1)
//p=p*r
tmp#=p.a
p.a=p.a*pt.a-p.b*pt.b
p.b=tmp#*pt.b+p.b*pt.a
next j
//f'(x) (=d) doesn't use the first value of cpoly() because the derivative of a constant
//is zero. So this iterated the above 1 more time just so the last value is accounted for.
f.a=f.a+p.a*cpoly(array count(cpoly())).a-p.b*cpoly(array count(cpoly())).b
f.b=f.b+p.b*cpoly(array count(cpoly())).a+p.a*cpoly(array count(cpoly())).b
//apply newton's method:
//rnew=r-f/d=r-f(r)/f'(r) (newton's method)
d2#=d.a*d.a+d.b*d.b
retval.a=pt.a-(f.a*d.a+f.b*d.b)/d2#
retval.b=pt.b-(d.a*f.b-f.a*d.b)/d2#
endfunction
Each point undergoes 1 iteration of newton's method per frame, and they quickly converge to the polynomial's four roots. The picture in my first snippet graphs the details of this process simply by measuring how long it took for the graph to converge.
The areas of nonconvergence are really precise, so pretty much all points will converge in the code.