Imitations of classical Julia set using Newton's method

[gamma]

I have been exploring the world of fractals derived from Newton's method for finding roots (solutions) to equations. The following page contains insights into some mathematics behind Newton's method.

http://www.chiark.greenend.org.uk/~sgtatham/newton/

If you open the page and step into the topic of "imitations" you will see what I mean.

Author mentions that we could try to generate Newton's fractal that looks like good old Julia set J(z^2+c) by solving a differential equation:

z_new = z -  f(z)/f'(z) = z^2 + c

We wonder, for which equation Newton's equation will produce z^2+c. Simon Tatham gives the solution to that equation which I don't understand. He treats this like the simplest form of differential equation which can be solved in a single step - integration of left and right side:

df/f = -dz / (z^2-z+c)

Here I am rather puzzled - if I may explain. Simon says that he solved the integral to the right as an integral of a rational function (http://www.answers.com/topic/partial-fractions-in-integration). Somehow he arrives at this equation:

f(z)= (z-a)^(1/(b-a)) * (z-b)^(1/(a-b))

a, b are roots of z^2-z+c=0.
 
(Please note that derivative of this function f can be simplified further for final calculation in computer using a perfectly well described method in the precious section of Simon's web site. That part is not a problem.)

I used Matlab software to check on this work. Matlab treats the initial differential equation as homogeneous differential equation and gives:

dsolve('z-f/Df=z^2-z+c','z')
 
ans =
 
C1*exp(-1/(c-1)^(1/2)*atan((z-1)/(c-1)^(1/2)))

Or, to follow Simon's story and compute integral of -dz/(z^2-z+c)

int(-1/(-z+z^2+c),z)
 
ans =
 
-2/(4*c-1)^(1/2)*atan((2*z-1)/(4*c-1)^(1/2))

Therefore I conclude I don't know how to solve it Simon's way and I need you to explain to me!

~~~~~~~~~~~~~~~~

If you're still following the story....

roots of z^2-z+c and some c are:

c=complex(0.2, 0.3)
p=[1 -1 c]
roots(p)

So far, that was the Matlab way.

~~~~~~~~~~~~~~~~

Now if you are really brave and enlightened like me of course, you will type into Matlab some other differential equation to see how it works in general:

dsolve('z-f/Df=z^5-z+c','z')
 
ans =
 
C1*exp(Int(-1/(-2*z+z^5+c),z))

which means Matlab left this integral for us to solve...

int(-1/(-2*z+z^5+c),z)
 
ans =
 
-sum(1/(-2+5*_R^4)*log(z-_R),_R = RootOf(-2*_Z+_Z^5+c))
 
which means that the output is now unreadable to me. Even though it seems like an error, "RootOf" is a legal way to mention roots of some polynomial. Mathematica also produces strange output for that integral: try here... http://integrals.wolfram.com/index.jsp

My second question is therefore, why are these programs for mathematics like this? why don't they work well? How do I read output if its good?

[gamma]

The interesting thing is that it works as expected, although I can not precisely say what is expected to be in detail, how good imitation can get. If formulas are just slightly more complicated a whole range of programs won't display good results or usable results or results that I could interpret in my humble understanding (ultrafractal, matlab, mathematica...). What is the course of action then? Why is the differential equation used here true?


Here is the code for ultrafractal (create .ufm file).



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
NewtonM {
;
;
;   Here we imitate Julia(z^2 + c) using
;   Newton's method...
;
;
;
;
;

init:

  complex zold = (0,0)
  complex i = (0,1)
  complex z = @start

loop:

     zold=z
     z = z - 1/((-0.5086-0.0747*i)/(z-(1.4623-0.1413*i)) + (0.5086+0.0747*i)/(z-(-0.4623+0.1413*i))) + #pixel

bailout:

     |z - zold| > @bailout && |z|<200

default:
  title = "Imitation of J(z^2+c)"
  maxiter = 2000
  periodicity = 0
  center = (0,0)
  magn = 1.2

  param start
    caption = "Start Value"
    default = (0,0)
  endparam

  param bailout
    caption = "Bailout"
    default = 0.000001
$IFDEF VER40
    exponential = true
$ENDIF
    hint = "Bailout value; smaller values will cause more \
            iterations to be done for each point."
  endparam


switch:

  type = "NewtonJ"

  s = #pixel
  bailout = @bailout

}

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
NewtonJ {
;
;
;   Here we imitate Julia(z^2 + c) using
;   Newton's method...
;
;
;
;
;
;

init:

  complex zold = (0,0)
  complex i = (0,1)
  complex z = #pixel

loop:

     zold=z

     z = z - 1/((-0.5086-0.0747*i)/(z-(1.4623-0.1413*i)) + (0.5086+0.0747*i)/(z-(-0.4623+0.1413*i))) + @s

bailout:

     |z - zold| > @bailout && |z|<200

default:
  title = "Imitation of J(z^2+c)"
  helpfile = "Uf*.chm"
  helptopic = "Html\formulas\standard\nova.html"
  maxiter = 2000
  periodicity = 0
  center = (0,0)
  magn = 0.75

  param s
    caption = "Julia Seed"
    default = (0,0)
    hint = "This is the Julia seed, Use the F7 Switch feature in \
            combination with the MandelNewton formula to find interesting \
            values."
  endparam

  param bailout
    caption = "Bailout"
    default = 0.000001
$IFDEF VER40
    exponential = true
$ENDIF
    hint = "Bailout value; smaller values will cause more \
            iterations to be done for each point."
  endparam


switch:

  type = "NewtonM"

  bailout = @bailout

}

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end

[gamma]

I received a reply that clarifies few things.
http://www.chiark.greenend.org.uk/~sgtatham/newton/

The math is good. There are no approximations in essence. With differential equation:
f=f(z)
df/f = dz / (z - z^2 - c)

inserted into mathematical software package such as Matlab I get:

dsolve('Df=f*(1/((z-z^2-c)))','z')
ans =
C1*exp(-2/(4*c-1)^(1/2)*atan((2*z-1)/(4*c-1)^(1/2)))

The fractal would be the same... But if we just factorize the denominator....

Example of factorizing 3rd degree polynomial:

c=complex(0.2, 0.3)
p=[1 -1 -c]
roots(p)

...and the roots are a and b,
 
dsolve('Df=f*(1/((z-a)*(z-b)))','z')
ans =
C1*(z-a)^(1/(a-b))*(z-b)^(-1/(a-b))

Drastic difference!

Using partial fractions...
http://www.mathsrevision.net/alevel/pages.php?page=94
http://cnx.org/content/m2111/latest/
... and following other instructions we can solve manually differential equation and get same form every time for most problems.

Partial fractions in Matlab:

define point c and polynomials p and q

c=complex(x,y)

p=[1]
q=[-1 1 -c]

[R,P,K]=residue(p,q)

It is now in the form
p/q = R1/(z-P1) + R2/(z-P2) + ...

~~~~~~~~~~~~~~~~~~~~~~~~~

This was conversion of one formula into function for which Newton's method will produce the fractal that we had with the first formula. There are differences however. Inevitably now we are studying *Newton's method* and its dynamics. For some functions it works bad - e.g. diverging from solutions. Sort of, initial formula is maybe wrapped around its solutions to some extent of fractal dimension... I'm watching Simpsons instead of this from now on.

[HPDZ]

Sorry I didn't read this until today...I could have saved you some headache with the partial fraction integrals. But I'm glad you were able to work it out. The version with Exp and Atan ought to reduce to the version obtained with partial fraction integration via some obscure trigonometric identities and complex number expansions.

We can continue the analysis a little further to try to visualize the structure of f(z) so we can understand the behavior of Newton's method as it tries to locate the zeros. We start with the function



We defined a and b as the roots of so we have



We can arbitrarily set to be the root with the + radical and the root with the - radical since is invariant under exchange of and .

For brevity define



and then insert all this into the equation for f and we get



As a complex-valued function of a complex argument involving a complex exponent, this is not a simple function to visualize. If we just look at what happens for , we at least can see it approaches 1 from below as and approaches 1 from above as . It of course has a root at (this is what Neton's method will be trying to find), and it has a singularity at .

The function can have a non-real value for some values of z, depending on what s is. And s can itself be imaginary if , which further complicates the structure of f.

I don't have access to my Mathematica system at the moment but I might try plotting some illuminating views of this this later.

[cKleinhuis]

hi there,
gamma seems to be a great self-talker ... respect, this is the way to enlightment!

and nice to see that you make use of the math functions !

but i have to admit i have not read it alll ... but nice to see active members here that like to share information !

 

[HPDZ]

and nice to see that you make use of the math functions !

Yes, I had not done much with Tex before but there was just no way to communicate in this case without it. It turned out to be easier than I thought, and quite powerful. I wish the baselines of the text were better aligned and the font sizes matched better....

Anyway: Here's a couple of simple quick plots of this thing for real values of Z.

The horizontal axis is Z and the vertical axis is f(Z) with f as described in the previous posts.

First with c=0. This gives s=1.



This has a zero at z=-0.5. Newton's method will converge to it for all positive starting values less than 1.5, where the singularity is. Most negative values will shoot off to the right and get zorched into z=+infinity land, but a small range very close to the zero at z=-0.5 will converge. You could work this out by finding those values of z that have a tangent line that intersects the z-axis at z<1.5.

The next graph is with c=-1. Now you can see there is a gap where f(z) takes on an imaginary component. I don't have time right now to make a better plot of this in 3D. You can see most negative values again get launched into the twilight zone vortex that exists for all real numbers to the right of the singularity. Only a tiny range of z values very close to the zero have any hope of convergence.



Well, I wish I could explore this more but no time right now.