Fractals/Mathematics/Newton method

< Fractals < Mathematics

Introduction

types

Animation of Newton's method in case of real-valued function with one root

Newton method can be used for finding successively better approximations to one root (zero) x of function f :

x:f(x)=0\,.

If one wants to find another root must apply the method again with other initial point

Newton method can be applied to :

a real-valued function ^[3]
a complex-valued function
a system of equations

How to find all roots ? ^[4]^[5]

systems of nonlinear equations

One may also use Newton's method to solve systems of ^[6]

2 (non-linear) equations
with 2 complex variables : $x_{1},x_{2}$ ^[7]

 ${\begin{cases}F_{1}(x_{1},x_{2})=0\\F_{2}(x_{1},x_{2})=0\end{cases}}$

It can be expressed more sinthetically using vector notation :

 ${\mathbf {F}}({\mathbf {x}})=0$

where ${\mathbf {x}}$ is a vector of 2 complex variables :

 ${\mathbf {x}}={\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}$

and ${\mathbf {F}}$ is a vector of functions ( or is a vector-valued function of the vector variable x ) :

  ${\mathbf {F}}={\begin{bmatrix}F_{1}\\F_{2}\end{bmatrix}}$

Solution can be find by iteration :^[8]

 ${\mathbf {x}}^{(k+1)}={\mathbf {x}}^{(k)}+{\mathbf {s}},\ k=0,1,2,\ldots .$

where s is an increment ( which is now a vector of 2 components ):

  ${\mathbf {s}}={\mathbf {x}}^{(k+1)}-{\mathbf {x}}^{(k)}$

Increment can be computed by :

 ${\mathbf {s}}=J^{-1}({\mathbf {x}}^{(k)}){\mathbf {F}}({\mathbf {x}}^{(k)})$

Here $J$ is Jacobian matrix^[9] of ${\mathbf {F}}$ at ${\mathbf {x}}^{(k)}$ :

  $J_{\mathbf {F}}({\mathbf {x}}^{(k)})=(F'({\mathbf {x}}))_{ij}={\dfrac {\partial F_{i}}{\partial x_{j}}}({\mathbf {x}})={\begin{bmatrix}{\dfrac {\partial F_{1}}{\partial x_{1}}}&{\dfrac {\partial F_{1}}{\partial x_{2}}}\\[1em]{\dfrac {\partial F_{2}}{\partial x_{1}}}&{\dfrac {\partial F_{2}}{\partial x_{2}}}\end{bmatrix}}$

and $J^{-1}$ is inverse of the Jacobian matrix $J$ .

In case of

a small number of equations use the inverse matrix.
a big number of equations : rather than actually computing the inverse of this matrix, one can save time by solving the system of linear equations $J_{F}(x_{n})(x_{n+1}-x_{n})=-F(x_{n})\,\!$ for the unknown x_n+1 − x_n.

Algorithm in pseudocode : ^[10]

start with an initial approximation ( guess) ${\mathbf {x}}^{k};k=0$
repeat until convergence ( iterate)
- solve: ${\mathbf {s}}=J^{-1}({\mathbf {x}}^{(k)}){\mathbf {F}}({\mathbf {x}}^{(k)})$
- compute : ${\mathbf {x}}^{(k+1)}={\mathbf {x}}^{(k)}+{\mathbf {s}}$
- $k=k+1$

Stop criteria :

absolute error : $|s|<1.e-12$ ^[11]

What is needed ?

function f ( a differentiable function )
it's derivative f'
starting point x0 ( which should be in the basin of root's attraction )
stopping rule ( criteria to stop iteration )^[12]

stopping rule

$|x_{n+1}-x_{n}|<\epsilon \;$

Pitfalls

method has a cycle and never converge (
- inflection point ^[13]
- local minimum
method is diverging not converging
- local minimum
multiple root ( multiplicity of root > 1 , for example : f(x) = f'(x) = 0). Slow approximation, derivative tends to zero, trouble in division step ( do not divide by zero ). Use modified Newton's method
slow convergence^[14]
- bad value of the initial approximation ( an initial guess far from the root )
- a function whose derivative vanishes near a root, or whose second derivative is unbounded near a root

Function has inflection point ( method has a cycle and never converge)
Local minimum

Pseudocode

The following is an example of using the Newton's Method to help find a root of a function f which has derivative fprime.

The initial guess will be $x_{0}=1$ and the function will be $f(x)=x^{2}-2$ so that $f'(x)=2x$ .

Each new iterative of Newton's method will be denoted by x1. We will check during the computation whether the denominator (yprime) becomes too small (smaller than epsilon), which would be the case if $f'(x_{n})\approx 0$ , since otherwise a large amount of error could be introduced.

%These choices depend on the problem being solved
x0 = 1                      %The initial value
f = @(x) x^2 - 2            %The function whose root we are trying to find
fprime = @(x) 2*x           %The derivative of f(x)
tolerance = 10^(-7)         %7 digit accuracy is desired
epsilon = 10^(-14)          %Don't want to divide by a number smaller than this

maxIterations = 20          %Don't allow the iterations to continue indefinitely
haveWeFoundSolution = false %Were we able to find the solution to the desired tolerance? not yet.

for i = 1 : maxIterations 

    y = f(x0)
    yprime = fprime(x0)

    if(abs(yprime) < epsilon)                         %Don't want to divide by too small of a number
        fprintf('WARNING: denominator is too small\n')
        break;                                        %Leave the loop
    end

    x1 = x0 - y/yprime                                %Do Newton's computation
    
    if(abs(x1 - x0)/abs(x1) < tolerance)              %If the result is within the desired tolerance
        haveWeFoundSolution = true
        break;                                        %Done, so leave the loop
    end

    x0 = x1                                           %Update x0 to start the process again                  
    
end

if (haveWeFoundSolution) % We found a solution within tolerance and maximum number of iterations
    fprintf('The root is: %f\n', x1);
else %If we weren't able to find a solution to within the desired tolerance
    fprintf('Warning: Not able to find solution to within the desired tolerance of %f\n', tolerance);
    fprintf('The last computed approximate root was %f\n', x1)
end

code

c ^[15]
matlab : newtonm function ^[16]
Maxima CAS
scilab ^[17]
Mathematica^[18]

applications of Newton's method in case of iterated complex quadratic polynomials

recurrence relations

Define the iterated function $f^{n}$ by:

${\begin{aligned}f^{0}(z,c)&=z\\f^{1}(z,c)&=z^{2}+c\\f^{m+1}(z,c)&=f(f^{m}(z,c),c)\end{aligned}}$

Pick arbitrary names for the iterated function and it's derivatives These derivatives can be computed by recurrence relations.

arbitrary names	Initial conditions	Recurrence steps
$A_{m}=f^{m}$	$A_{0}=z$	$A_{m+1}=A_{m}^{2}+c$
$B_{m}={\dfrac {\partial }{\partial z}}f^{m}$	$B_{0}=1$	$B_{m+1}=2A_{m}B_{m}$
$C_{m}={\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}f^{m}$	$C_{0}=0$	$C_{m+1}=2(B_{m}^{2}+A_{m}C_{m})$
$D_{m}={\dfrac {\partial }{\partial c}}f^{m}$	$D_{0}=0$	$D_{m+1}=2A_{m}D_{m}+1$
$E_{m}={\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}f^{m}$	$E_{0}=0$	$E_{m+1}=2(A_{m}E_{m}+B_{m}D_{m})$

Derivation of recurrence relations :

$A_{0}\quad {\xrightarrow {definition\ of\ A}}\ f^{0}(z,c)\quad {\xrightarrow {definition\ of\ f}}\ \quad z$

$B_{0}\quad {\xrightarrow {definition\ of\ B}}\ {\dfrac {\partial }{\partial z}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}\ {\dfrac {\partial }{\partial z}}z{\xrightarrow {derivative}}\ 1$

$C_{0}\quad {\xrightarrow {definition\ of\ C}}\ {\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}\ {\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}z{\xrightarrow {derivative}}{\dfrac {\partial }{\partial z}}1{\xrightarrow {derivative}}0$

$D_{0}\quad {\xrightarrow {definition\ of\ D}}{\dfrac {\partial }{\partial c}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}{\dfrac {\partial }{\partial c}}z{\xrightarrow {derivative}}0$

$E_{0}\quad {\xrightarrow {definition\ of\ E}}{\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}f^{0}(z,c)\quad {\xrightarrow {definitionoff}}{\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}z\quad {\xrightarrow {derivative}}{\dfrac {\partial }{\partial c}}1\quad {\xrightarrow {derivative}}0$

$A_{m+1}\quad {\xrightarrow {definitionofA}}f^{m+1}(z,c)\quad {\xrightarrow {definition\ of\ f}}f(f^{m}(z,c),c)\quad {\xrightarrow {definition\ of\ A}}f(A_{m},c)\quad {\xrightarrow {definition\ of\ f}}A_{m}^{2}+c$

$B_{m+1}\quad {\xrightarrow {definition\ of\ B}}{\dfrac {\partial }{\partial z}}A_{m+1}\quad {\xrightarrow {definition\ of\ A}}{\dfrac {\partial }{\partial z}}(A_{m}^{2}+c)\quad {\xrightarrow {distributivity}}{\dfrac {\partial }{\partial z}}A_{m}^{2}+{\dfrac {\partial }{\partial z}}c\quad {\xrightarrow {constant\ derivative}}{\dfrac {\partial }{\partial z}}A_{m}^{2}+0\quad {\xrightarrow {zero}}{\dfrac {\partial }{\partial z}}A_{m}^{2}\quad {\xrightarrow {chain\ rule}}2A_{m}({\dfrac {\partial }{\partial z}}A_{m})\quad {\xrightarrow {definition\ of\ B}}2A_{m}B_{m}$

$C_{m+1}\quad {\xrightarrow {definition\ of\ C}}{\dfrac {\partial }{\partial z}}B_{m+1}\quad {\xrightarrow {definition\ of\ B}}{\dfrac {\partial }{\partial z}}(2A_{m}B_{m})\quad {\xrightarrow {linearity}}2{\dfrac {\partial }{\partial z}}(A_{m}B_{m})\quad {\xrightarrow {product\ rule}}2(({\dfrac {\partial }{\partial z}}A_{m})B_{m}+A_{m}({\dfrac {\partial }{\partial z}}B_{m}))\quad {\xrightarrow {definition\ of\ B}}2(B_{m}B_{m}+A_{m}({\dfrac {\partial }{\partial z}}B_{m}))\quad {\xrightarrow {algebra}}2(B_{m}^{2}+A_{m}({\dfrac {\partial }{\partial z}}B_{m}))\quad {\xrightarrow {definition\ of\ C}}2(B_{m}^{2}+A_{m}C_{m})$

$D_{m+1}\quad {\xrightarrow {definition\ of\ D}}{\dfrac {\partial }{\partial c}}A_{m+1}\quad {\xrightarrow {definition\ of\ A}}{\dfrac {\partial }{\partial c}}(A_{m}^{2}+c)\quad {\xrightarrow {distributivity}}{\dfrac {\partial }{\partial c}}A_{m}^{2}+{\dfrac {\partial }{\partial c}}c\quad {\xrightarrow {derivative}}{\dfrac {\partial }{\partial c}}A_{m}^{2}+1\quad {\xrightarrow {chain\ rule}}2A_{m}({\dfrac {\partial }{\partial c}}A_{m})+1\quad {\xrightarrow {definition\ of\ D}}2A_{m}D_{m}+1$

$E_{m+1}\quad {\xrightarrow {definition\ of\ E}}{\dfrac {\partial }{\partial c}}B_{m+1}\quad {\xrightarrow {definition\ of\ B}}{\dfrac {\partial }{\partial c}}(2A_{m}B_{m})\quad {\xrightarrow {linearity}}2{\dfrac {\partial }{\partial c}}(A_{m}B_{m})\quad {\xrightarrow {product\ rule}}2(A_{m}({\dfrac {\partial }{\partial c}}B_{m})+({\dfrac {\partial }{\partial c}}A_{m})B_{m})\quad {\xrightarrow {definitionof\ E}}2(A_{m}E_{m}+({\dfrac {\partial }{\partial c}}A_{m})B_{m})\quad {\xrightarrow {definition\ of\ D}}2(A_{m}E_{m}+D_{m}B_{m})$

center

To use it in a GUI program :

click on the menu item : nucleus
using a mouse mark rectangle region around center of hyperbolic component

A center ( or nucleus ) $c=n$ of period $p$ satisfies :

f^{p}(0,n)=0

Applying Newton's method in one variable:

$n_{m+1}=N(n_{m})=n_{m}-{\frac {f^{p}(0,n_{m})}{{\dfrac {\partial }{\partial c}}f^{p}(0,n_{m})}}=n_{m}-{\frac {A_{p}(0,n_{m})}{D_{p}(0,n_{m})}}$

where initial estimate $n_{0}$ is arbitrary choosen ( abs(n) <2.0 )

stoppping rule : $A_{p}=0\Leftrightarrow n_{m+1}=n_{m}$

arbitrary names	Initial conditions	Recurrence steps
$A_{m}=f^{m}$	$A_{0}=z=0$	$A_{m+1}=A_{m}^{2}+c$
$D_{m}={\dfrac {\partial }{\partial c}}f^{m}$	$D_{0}=0$	$D_{m+1}=2A_{m}D_{m}+1$

/*


code from  unnamed c program ( book ) 
http://code.mathr.co.uk/book
see mandelbrot_nucleus.c

by Claude Heiland-Allen




COMPILE :

 gcc -std=c99 -Wall -Wextra -pedantic -O3 -ggdb  m.c -lm -lmpfr

usage: 

./a.out bits cx cy period maxiters
    
output :  space separated complex nucleus on stdout
    
example 

./a.out 53 -1.75 0 3 100
./a.out 53 -1.75 0 30 100
./a.out 53  0.471 -0.3541 14 100
./a.out 53 -0.12 -0.74 3 100
*/



#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h> // arbitrary precision
#include <mpfr.h>




extern int mandelbrot_nucleus(mpfr_t cx, mpfr_t cy, const mpfr_t c0x, const mpfr_t c0y, int period, int maxiters) {
  int retval = 0;
  mpfr_t zx, zy, dx, dy, s, t, u, v;
  mpfr_inits2(mpfr_get_prec(c0x), zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);
  mpfr_set(cx, c0x, GMP_RNDN);
  mpfr_set(cy, c0y, GMP_RNDN);

  for (int i = 0; i < maxiters; ++i) {
    // z = 0
    mpfr_set_ui(zx, 0, GMP_RNDN);
    mpfr_set_ui(zy, 0, GMP_RNDN);
    // d = 0
    mpfr_set_ui(dx, 0, GMP_RNDN);
    mpfr_set_ui(dy, 0, GMP_RNDN);
    for (int p = 0; p < period; ++p) {
      // d = 2 * z * d + 1;
      mpfr_mul(u, zx, dx, GMP_RNDN);
      mpfr_mul(v, zy, dy, GMP_RNDN);
      mpfr_sub(u, u, v, GMP_RNDN);
      mpfr_mul_2ui(u, u, 1, GMP_RNDN);
      mpfr_mul(dx, dx, zy, GMP_RNDN);
      mpfr_mul(dy, dy, zx, GMP_RNDN);
      mpfr_add(dy, dx, dy, GMP_RNDN);
      mpfr_mul_2ui(dy, dy, 1, GMP_RNDN);
      mpfr_add_ui(dx, u, 1, GMP_RNDN);
      // z = z^2 + c;
      mpfr_sqr(u, zx, GMP_RNDN);
      mpfr_sqr(v, zy, GMP_RNDN);
      mpfr_mul(zy, zx, zy, GMP_RNDN);
      mpfr_sub(zx, u, v, GMP_RNDN);
      mpfr_mul_2ui(zy, zy, 1, GMP_RNDN);
      mpfr_add(zx, zx, cx, GMP_RNDN);
      mpfr_add(zy, zy, cy, GMP_RNDN);
    }
    // check d == 0
    if (mpfr_zero_p(dx) && mpfr_zero_p(dy)) {
      retval = 1;
      goto done;
    }

    
    // st = c - z / d
    mpfr_sqr(u, dx, GMP_RNDN);
    mpfr_sqr(v, dy, GMP_RNDN);
    mpfr_add(u, u, v, GMP_RNDN);
    mpfr_mul(s, zx, dx, GMP_RNDN);
    mpfr_mul(t, zy, dy, GMP_RNDN);
    mpfr_add(v, s, t, GMP_RNDN);
    mpfr_div(v, v, u, GMP_RNDN);
    mpfr_mul(s, zy, dx, GMP_RNDN);
    mpfr_mul(t, zx, dy, GMP_RNDN);
    mpfr_sub(zy, s, t, GMP_RNDN);
    mpfr_div(zy, zy, u, GMP_RNDN);
    mpfr_sub(s, cx, v, GMP_RNDN);
    mpfr_sub(t, cy, zy, GMP_RNDN);
    // uv = st - c
    mpfr_sub(u, s, cx, GMP_RNDN);
    mpfr_sub(v, t, cy, GMP_RNDN);
    // c = st
    mpfr_set(cx, s, GMP_RNDN);
    mpfr_set(cy, t, GMP_RNDN);
    // check uv = 0
    if (mpfr_zero_p(u) && mpfr_zero_p(v)) {
      retval = 2;
      goto done;
    }
  } // for (int i = 0; i < maxiters; ++i)


 done: mpfr_clears(zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);


 return retval;
}



void DescribeStop(int stop)
{
  switch( stop )
  {
    case 0:
    printf(" method stopped because i = maxiters\n");
    break;
    
    case 1:
    printf(" method stopped because derivative == 0\n");
    break;
    
    //...
    case 2:
    printf(" method stopped because uv = 0\n");
    break;
    
   }
}



void usage(const char *progname) {
  fprintf(stderr,
    "program finds one center ( nucleus) of hyperbolic component of Mandelbrot set using Newton method"
    "usage: %s bits cx cy period maxiters\n"
    "outputs space separated complex nucleus on stdout\n"
    "example %s 53 -1.75 0 3 100\n",
    progname, progname);
}




int main(int argc, char **argv) {
  
  // check the input 
  if (argc != 6) { usage(argv[0]); return 1; }
  
  // read the values 
  int bits = atoi(argv[1]);
  mpfr_t cx, cy, c0x, c0y;
  mpfr_inits2(bits, cx, cy, c0x, c0y, (mpfr_ptr) 0);
  mpfr_set_str(c0x, argv[2], 10, GMP_RNDN);
  mpfr_set_str(c0y, argv[3], 10, GMP_RNDN);
  int period = atoi(argv[4]);
  int maxiters = atoi(argv[5]);
  int stop; 


  //
  stop = mandelbrot_nucleus(cx, cy, c0x, c0y, period, maxiters);

   
  //
  printf(" nucleus ( center) of component with period = %s near c = %s ; %s is : \n ", argv[4], argv[2], argv[3]);
  mpfr_out_str(0, 10, 0, cx, GMP_RNDN);
  putchar(' ');
  mpfr_out_str(0, 10, 0, cy, GMP_RNDN);
  putchar('\n');
  //
  DescribeStop(stop) ;

  // clear memeory 
  mpfr_clears(cx, cy, c0x, c0y, (mpfr_ptr) 0);
  
  // 
  return 0;
}

boundary

The boundary of the component with center $n$ of period $p$ can be parameterized by internal angles.

The boundary point $c=b$ at angle $t$ (measured in turns) satisfies system of 2 equations :

${\begin{cases}f^{p}(w,b)-w&=0\\{\dfrac {\partial }{\partial z}}f^{p}(w,b)-e^{2\pi it}&=0\end{cases}}$

Defining a function of two complex variables:

${\begin{aligned}g{\begin{pmatrix}z\\c\end{pmatrix}}&={\begin{pmatrix}f^{p}(z,c)-z\\{\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it}\end{pmatrix}}\end{aligned}}$

and applying Newton's method in two variables:

${\begin{aligned}{\boldsymbol {v}}_{0}={\begin{pmatrix}n\\n\end{pmatrix}}\\J_{g}({\boldsymbol {v}}_{m})({\boldsymbol {v}}_{m+1}-{\boldsymbol {v}}_{m})=-g({\boldsymbol {v}}_{m})\end{aligned}}$

where

${\begin{aligned}J_{g}={\begin{pmatrix}{\dfrac {\partial }{\partial z}}(f^{p}(z,c)-z)&{\dfrac {\partial }{\partial c}}(f^{p}(z,c)-z)\\{\dfrac {\partial }{\partial z}}({\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it})&{\dfrac {\partial }{\partial c}}({\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it})\end{pmatrix}}\end{aligned}}$

This can be expressed using the recurrence relations as:

${\begin{aligned}{\begin{pmatrix}w_{0}\\b_{0}\end{pmatrix}}&={\begin{pmatrix}n\\n\end{pmatrix}}\\{\begin{pmatrix}B_{p}-1&D_{p}\\C_{p}&E_{p}\end{pmatrix}}{\begin{pmatrix}w_{m+1}-w_{m}\\b_{m+1}-b_{m}\end{pmatrix}}&=-{\begin{pmatrix}A_{p}-w_{m}\\B_{p}-e^{2\pi it}\end{pmatrix}}\end{aligned}}$

where $\{A,B,C,D,E\}_{p}$ are evaluated at $(w_{m},b_{m})$ .

Example

"The process is for numerically finding one boundary point on one component at a time - you can't solve for the whole boundary of all the components at once. So pick and fix one particular nucleus n and one particular internal angle t " (Claude Heiland-Allen ^[19])

Lets try find point of boundary ( bond point )

$c=b$

of hyperbolic component of Mandelbrot set for :

period p=3
angle t=0
center ( nucleus ) c = n = -0.12256116687665+0.74486176661974i

There is only one such point.

Initial estimates (first guess) will be center ( nucleus ) $c=n$ of this components. This center ( complex number) will be saved in a vector of complex numbers containing two copies of that nucleus:

${\begin{aligned}{\boldsymbol {v}}_{0}={\begin{pmatrix}n\\n\end{pmatrix}}={\begin{pmatrix}-0.12256116687665+0.74486176661974i\\-0.12256116687665+0.74486176661974i\end{pmatrix}}\end{aligned}}$

The boundary point at angle $t$ (measured in turns) satisfies system of 2 equations :

${\begin{cases}f^{3}(z,c)-z&=0\\{\dfrac {\partial }{\partial z}}f^{3}(z,c)-e^{2\pi i0}&=0\end{cases}}$

so :

${\begin{cases}z^{8}+4cz^{6}+6c^{2}z^{4}+2cz^{4}+4c^{3}z^{2}+4c^{2}z^{2}-z+c^{4}+2c^{3}+c^{2}+c=0\\8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1=0\end{cases}}$

Then function g

${\begin{aligned}g({\boldsymbol {v}})={\begin{cases}z^{8}+4cz^{6}+6c^{2}z^{4}+2cz^{4}+4c^{3}z^{2}+4c^{2}z^{2}-z+c^{4}+2c^{3}+c^{2}+c\\8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1\end{cases}}\end{aligned}}$

Jacobian matrix $J$ is

$J_{g}({\boldsymbol {v}})={\begin{pmatrix}8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1&4z^{6}+12cz^{4}+2z^{4}+12c^{2}z^{2}+8cz^{2}+4c^{3}+6c^{2}+2c+1\\56z^{6}+120cz^{4}+72c^{2}z^{2}+24cz^{2}+8c^{3}+8c^{2}&24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz\end{pmatrix}}$

$J_{g}^{-1}({\boldsymbol {v}})={\begin{pmatrix}{\frac {24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz}{d}}&{\frac {-4*z^{6}-12*c*z^{4}-2*z^{4}-12*c^{2}*z^{2}-8*c*z^{2}-4*c^{3}-6*c^{2}-2*c-1}{d}}\\{\frac {-56*z^{6}-120*c*z^{4}-72*c^{2}*z^{2}-24*c*z^{2}-8*c^{3}-8*c^{2}}{d}}&{\frac {8*z^{7}+24*c*z^{5}+24*c^{2}*z^{3}+8*c*z^{3}+8*c^{3}*z+8*c^{2}*z-1}{d}}\end{pmatrix}}$

where denominator d is :

$d=(24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz)(8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1)+(-56z^{6}-120cz^{4}-72c^{2}z^{2}-24cz^{2}-8c^{3}-8c^{2})(4z^{6}+12cz^{4}+2z^{4}+12c^{2}z^{2}+8cz^{2}+4c^{3}+6c^{2}+2c+1)$

Then find better aproximation of point c=b using iteration :

${\boldsymbol {v}}_{m+1}={\frac {{\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m})}{J_{g}({\boldsymbol {v}}_{m})}}=({\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m}))J_{g}^{-1}({\boldsymbol {v}}_{m})$

using ${\boldsymbol {v}}_{0}$ as an starting point :

${\boldsymbol {v}}_{0}$

${\boldsymbol {v}}_{1}=({\boldsymbol {v}}_{0}-g({\boldsymbol {v}}_{0}))J_{g}^{-1}({\boldsymbol {v}}_{0})$

${\boldsymbol {v}}_{2}=({\boldsymbol {v}}_{1}-g({\boldsymbol {v}}_{1}))J_{g}^{-1}({\boldsymbol {v}}_{1})$

$...$

${\boldsymbol {v}}_{m+1}=({\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m}))J_{g}^{-1}({\boldsymbol {v}}_{m})$

Maxima CAS code

Boundaries of Components of Mandelbrot set by Newton method. Image and Maxima CAS code

C++ code

/*
cpp code by Claude Heiland-Allen 
http://mathr.co.uk/blog/
licensed GPLv3+


g++ -std=c++98 -Wall -pedantic -Wextra -O3 -o bonds bonds.cc

./bonds period nucleus_re nucleus_im internal_angle
./bonds 1 0 0 0
./bonds 1 0 0 0.5
./bonds 1 0 0 0.333333333333333
./bonds 2 -1 0 0
./bonds 2 -1 0 0.5
./bonds 2 -1 0 0.333333333333333
./bonds 3 -0.12256116687665 0.74486176661974 0
./bonds 3 -0.12256116687665 0.74486176661974 0.5
./bonds 3 -0.12256116687665 0.74486176661974 0.333333333333333

for b in $(seq -w 0 999)
do
  ./bonds 1 0 0 0.$b 2>/dev/null
done > period-1.txt
gnuplot
plot "./period-1.txt" with lines


-------------
for b in $(seq -w 0 999)
do
  ./bonds 3 -0.12256116687665 0.74486176661974 0.$b 2>/dev/null
done > period-3a.txt

gnuplot
plot "./period-3a.txt" with lines
*/
#include <complex>
#include <stdio.h>
#include <stdlib.h>

typedef std::complex<double> complex;

const double pi = 3.141592653589793;
const double eps = 1.0e-12;

struct param {
  int period;
  complex nucleus, bond;
};

struct recurrence {
  complex A, B, C, D, E;
};

void recurrence_init(recurrence &r, const complex &z) {
  r.A = z;
  r.B = 1;
  r.C = 0;
  r.D = 0;
  r.E = 0;
}

/*
void recurrence_step(recurrence &rr, const recurrence &r, const complex &c) {
  rr.A = r.A * r.A + c;
  rr.B = 2.0 * r.A * r.B;
  rr.C = 2.0 * (r.B * r.B + r.A * r.C);
  rr.D = 2.0 * r.A * r.D + 1.0;
  rr.E = 2.0 * (r.A * r.E + r.B * r.D);
}
*/

void recurrence_step_inplace(recurrence &r, const complex &c) {
  // this works because no component of r is read after it is written
  r.E = 2.0 * (r.A * r.E + r.B * r.D);
  r.D = 2.0 * r.A * r.D + 1.0;
  r.C = 2.0 * (r.B * r.B + r.A * r.C);
  r.B = 2.0 * r.A * r.B;
  r.A = r.A * r.A + c;
}

struct matrix {
  complex a, b, c, d;
};

struct vector {
  complex u, v;
};

vector operator+(const vector &u, const vector &v) {
  vector vv = { u.u + v.u, u.v + v.v };
  return vv;
}

vector operator*(const matrix &m, const vector &v) {
  vector vv = { m.a * v.u + m.b * v.v, m.c * v.u + m.d * v.v };
  return vv;
}

matrix operator*(const complex &s, const matrix &m) {
  matrix mm = { s * m.a, s * m.b, s * m.c, s * m.d };
  return mm;
}

complex det(const matrix &m) {
  return m.a * m.d - m.b * m.c;
}

matrix inv(const matrix &m) {
  matrix mm = { m.d, -m.b, -m.c, m.a };
  return (1.0/det(m)) * mm;
}

// newton stores <z, c> into vector as <u, v>

vector newton_init(const param &h) {
  vector vv = { h.nucleus, h.nucleus };
  return vv;
}

void print_equation(const matrix &m, const vector &v, const vector &k) {
  fprintf(stderr, "/  % 20.16f%+20.16fi  % 20.16f%+20.16fi  \\  /  z  -  % 20.16f%+20.16fi  \\  =  /  % 20.16f%+20.16fi  \\\n", real(m.a), imag(m.a), real(m.b), imag(m.b), real(v.u), imag(v.u), real(k.u), imag(k.u));
  fprintf(stderr, "\\  % 20.16f%+20.16fi  % 20.16f%+20.16fi  /  \\  c  -  % 20.16f%+20.16fi  /     \\  % 20.16f%+20.16fi  /\n", real(m.c), imag(m.c), real(m.d), imag(m.d), real(v.v), imag(v.v), real(k.v), imag(k.v));
}

vector newton_step(const vector &v, const param &h) {
  recurrence r;
  recurrence_init(r, v.u);
  for (int i = 0; i < h.period; ++i) {
    recurrence_step_inplace(r, v.v);
  }
  // matrix equation: J * (vv - v) = -g
  // solved by: vv = inv(J) * -g + v
  // note: the C++ variable g has the value of semantic variable -g
  matrix J = { r.B - 1.0, r.D, r.C, r.E };
  vector g = { v.u - r.A, h.bond - r.B };
  print_equation(J, v, g);
  return inv(J) * g + v;
}

int main(int argc, char **argv) {
  if (argc < 5) {
    fprintf(stderr, "usage: %s period nucleus_re nucleus_im internal_angle\n", argv[0]);
    return 1;
  }
  param h;
  h.period  = atoi(argv[1]);
  h.nucleus = complex(atof(argv[2]), atof(argv[3]));
  double angle = atof(argv[4]);
  if (angle == 0 && h.period != 1) {
    fprintf(stderr, "warning: possibly wrong results due to convergence into parent!\n");
  }
  h.bond = std::polar(1.0, 2.0 * pi * angle);
  fprintf(stderr, "initial parameters:\n");
  fprintf(stderr, "period = %d\n", h.period);
  fprintf(stderr, "nucleus = % 20.16f%+20.16fi\n", real(h.nucleus), imag(h.nucleus));
  fprintf(stderr, "bond = % 20.16f%+20.16fi\n", real(h.bond), imag(h.bond));
  vector v = newton_init(h);
  fprintf(stderr, "z =% 20.16f%+20.16fi\n", real(v.u), imag(v.u));
  fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
  fprintf(stderr, "\n");
  fprintf(stderr, "newton iteration:\n");
  vector vv;
  do {
    vv = v;
    v = newton_step(vv, h);
    fprintf(stderr, "z =% 20.16f%+20.16fi\n", real(v.u), imag(v.u));
    fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
    fprintf(stderr, "\n");
  } while (abs(v.u - vv.u) + abs(v.v - vv.v) > eps);
  fprintf(stderr, "bond location:\n");
  fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
  // suitable for gnuplot
  printf("%.16f %.16f\n", real(v.v), imag(v.v));
  return 0;
}

size

Newton's method (providing the initial estimate is sufficiently good and the function is well-behaved) provides successively more accurate approximations. How to tell when the approximation is good enough? One way might be to compare the center with points on its boundary, and continue to increase precision until this distance is accurate to enough bits.

Algorithm:

given a center location estimate $n_{m}$ accurate to $P$ bits of precision
compute a boundary location estimate $b_{m}$ accurate to the same precision, using center $n_{m}$ as starting value
compute $|n_{m}-b_{m}|$ to find an estimate of the size of the component
1. if it isn't zero, and if it is accurate enough (to a few 10s of bits, perhaps) then finish
2. otherwise increase precision, refine center estimate to new precision, try again from the top

Measuring effective accuracy of the distance between two very close points might be done by comparing floating point exponents with floating point mantissa precision.

internal coordinate

https://mathr.co.uk/blog/2013-04-01_interior_coordinates_in_the_mandelbrot_set.html

It is possible to map hyperbolic component to closed unit disk centered at origin :

${\bar {D}}_{1}=\{b:\left|b\right|\leq 1\}$ .

This relation is described by system of 2 equations :

 ${\begin{cases}F^{p}(z,c)=z\\{\frac {\partial }{\partial z}}F^{p}(z,c)=b\end{cases}}$

where

p is the period of the target hyperbolic component on parameter plane
$\theta$ the desired internal angle of points c and b
$r \le 1$ is an internal radius
$b=r*e^{2\pi i\theta }$ is a point of unit circle = internal coordinate
c is a point of parameter plane
z is a point of dynamic plane
$(F^{0}(z,c)=z$ and $F^{q+1}(z,c)=F^{q}(F(z,c)^{2}+c)$

The algorithm by Claude Heiland-Allen for finding internal coordinate b from c :

choose c
check c ( bailout test on dynamical plane). When c is outside the Mandelbrot set, give up now ( or compute external coordinate );
Start with period one : $p=1$
while $p<p_{max}$
- Find $z_{0}$ such that $F^{p}(z_{0},c)=z_{0}$ using Newton's method in one complex variable;
- Find b by evaluating ${\frac {\partial }{\partial z}}F^{p}(z,c)$ at $z_{0}$
- If $\left|b\right|\leq 1$ then return b
- otherwise continue with the next p : $p=p+1$

To solve such system of equations one can use Newton method.

internal ray

Green internal and red external rays

/*
http://code.mathr.co.uk/mandelbrot-graphics/blob/HEAD:/c/bin/m-subwake-diagram-a.c
by Claude Heiland-Allen
*/

#include <mandelbrot-graphics.h>
#include <mandelbrot-numerics.h>
#include <mandelbrot-symbolics.h>
#include <cairo.h>

const double twopi = 6.283185307179586;


void draw_internal_ray(m_image *image, m_d_transform *transform, int period, double _Complex nucleus, const char *angle, double pt, m_pixel_t colour) {
  int steps = 128;
  mpq_t theta;
  mpq_init(theta);
  mpq_set_str(theta, angle, 10);
  mpq_canonicalize(theta);
  double a = twopi * mpq_get_d(theta);
  mpq_clear(theta);
  double _Complex interior = cos(a) + I * sin(a);

  double _Complex cl = 0, cl2 = 0;
  double _Complex c = nucleus;
  double _Complex z = c;
  cairo_surface_t *surface = m_image_surface(image);
  cairo_t *cr = cairo_create(surface);
  cairo_set_source_rgba(cr, m_pixel_red(colour), m_pixel_green(colour), m_pixel_blue(colour), m_pixel_alpha(colour));
  for (int i = 0; i < steps; ++i) {
    if (2 * i == steps) {
      cl = c;
    }
    if (2 * i == steps + 2) {
      cl2 = c;
    }
    double radius = (i + 0.5) / steps;
    m_d_interior(&z, &c, z, c, radius * interior, period, 64);
    double _Complex pc = c;
    double _Complex pdc = 1;
    m_d_transform_reverse(transform, &pc, &pdc);
    if (i == 0) {
      cairo_move_to(cr, creal(pc), cimag(pc));
    } else {
      cairo_line_to(cr, creal(pc), cimag(pc));
    }
  }
  cairo_stroke(cr);
  if (a != 0) {
    double t = carg(cl2 - cl);
    cairo_save(cr);
    double _Complex dcl = 1;
    m_d_transform_reverse(transform, &cl, &dcl);
    cairo_translate(cr, creal(cl), cimag(cl));
    cairo_rotate(cr, -t);
    cairo_translate(cr, 0, -pt/3);
    cairo_select_font_face(cr, "LMSans10", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL);
    cairo_set_font_size(cr, pt);
    cairo_text_path(cr, angle);
    cairo_fill(cr);
    cairo_restore(cr);
  }
  cairo_destroy(cr);
}

external rays

Dynamic ray

//  qmnplane.cpp by Wolf Jung (C) 2007-2014.
// part of Mandel 5.11 http://mndynamics.com/indexp.html
// the GNU General   Public License

//Time ~ nmax^2 , therefore limited  nmax .
int QmnPlane::newtonRay(int signtype, qulonglong N, qulonglong D,
   double &a, double &b, int q, QColor color, int mode) //5, white, 1
{  uint n; int j, i, k, i0, k0; mndAngle t; t.setAngle(N, D, j);
   double logR = 14.0, x, y, u;
   u = exp(0.5*logR); y = t.radians();
   x = u*cos(y); y = u*sin(y);
   if (pointToPos(x, y, i0, k0) > 1) i0 = -10;
   stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
   for (n = 1; n <= (nmax > 5000u ? 5000u : nmax + 2); n++)
   {  t.twice();
      for (j = 1; j <= q; j++)
      {  if (mode & 4 ? tricornNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians())
           : rayNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians()) )
         { n = (n <= 20 ? 65020u : 65010u); break; }
         if (pointToPos(x, y, i, k) > 1) i = -10;
         if (i0 > -10)
         {  if (i > -10) { if (!(mode & 8)) p->drawLine(i0, k0, i, k); }
            else { n = 65020u; break; }
         }
         i0 = i; k0 = k;
      }
   }
   //if rayNewton fails after > 20 iterations, endpoint may be accepted
   delete p; update(); if (n >= 65020u) return 2;
   if (mode & 2) { a = x; b = y; }
   if (n >= 65010u) return 1; else return 0;
} //newtonRay

int QmnPlane::rayNewton(int signtype, uint n, double a, double b,
   double &x, double &y, double rlog, double ilog)
{  uint k, l; double fx, fy, px, py, u, v = 0.0, d = 1.0 + x*x + y*y, t0, t1;
   for (k = 1; k <= 60; k++)
   {  if (signtype > 0) { a = x; b = y; }
      fx = cos(ilog); fy = sin(ilog);
      t0 = exp(rlog)*fx - 0.5*exp(-rlog)*(a*fx + b*fy);
      t1 = exp(rlog)*fy + 0.5*exp(-rlog)*(a*fy - b*fx);
      fx = x; fy = y; px = 1.0; py = 0.0;
      for (l = 1; l <= n; l++)
      {  u = 2.0*(fx*px - fy*py); py = 2.0*(fx*py + fy*px);
         px = u; if (signtype > 0) px++;
         u = fx*fx; v = fy*fy; fy = 2.0*fx*fy + b; fx = u - v + a;
         u += v; v = px*px + py*py; if (u + v > 1.0e100) return 1;
      }
      fx -= t0; fy -= t1; if (v < 1.0e-50) return 2;
      u = (fx*px + fy*py)/v; v = (fx*py - fy*px)/v;
      px = u*u + v*v; if (px > 9.0*d) return -1;
      x -= u; y += v; d = px; if (px < 1.0e-28 && k >= 5) break;
   }
   return 0;
} //rayNewton

Parameter ray

External parameter rays

double externalAngle(...) {
  ... 
  return (std::atan2(cky,ckx));
}

"This gets you the angle in only double-precision, but using double precision floating point throughout it's possible to get the external angle in much higher precision - the trick is to collect bits from the binary representation of the angle as you cross each dwell band - whether the final iterate that escaped has a positive or negative imaginary part determines if the bit is 0 or 1 respectively, see binary decomposition colouring http://www.mrob.com/pub/muency/binarydecomposition.html .

This is orthogonal to using Newton's method - it should also work with your potential field particle simulation technique as long as you never skip more than one dwell band at a time.

The idea behind Newton's method is to gradually reduce your target potential keeping the angle fixed to find a target final-z value, then solve for c in f^n(0, c) = z-final using the current c-value as initial guess for Newton's method iterations. The amount you reduce the potential by is controlled by a parameter called "sharpness", essentially the number of steps to take within each dwell band - usually I use between 4 and 8. This means the number of steps taken is known in advance (sharpness * initial dwell) and usually you don't need many Newton's method iterations to converge at each step (around 2 - 5 was typical in some tests I did a while ago).

You can only decrease potential so much for a fixed n with a fixed escape radius, so eventaully you decrement n, which increases |z-final| and gives you a choice of two angles (the pre-images under angle doubling) for arg z-final, you pick the one which is closer to the n'th iterate of your current c value, and which one you pick also tells you another bit of the external angle).

The algorithm is described here, for the other direction (starting from an external angle and tracing inwards to the Mandelbrot set - it's simpler because there is only one image under angle-doubling instead of two possibilities for the inverse): http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf

I have a few implementations, the currently maintained ones are:

external ray tracing inwards (double precision floating point, libgmp arbitrary precision rational for external angle) https://gitorious.org/mandelbrot/mandelbrot-numerics/source/3d55cfc99cc97decb0ba0d9c2fb271a504b8e504:c/lib/m_d_exray_in.c

external ray tracing outwards (double precision floating point, libgmp arbitrary precision rational for external angle) https://gitorious.org/mandelbrot/mandelbrot-numerics/source/3d55cfc99cc97decb0ba0d9c2fb271a504b8e504:c/lib/m_d_exray_out.c

external ray tracing inwards (libmpfr arbitrary precision floating point, libgmp arbitrary precision rational for external angle) https://gitorious.org/mandelbrot/mandelbrot-numerics/source/3d55cfc99cc97decb0ba0d9c2fb271a504b8e504:c/lib/m_r_exray_in.c " Claude Heiland-Allen on Fractal Forum^[20]

Compute one point of the ray

What one needs to start :

external angle $\theta$ of the ray on wants to draw. Angle is usually given in turns
function P ( which approximates Boettcher mapping )
it's derivative P'
starting point c0
stopping rule ( criteria to stop iteration )

arbitrary names	Initial conditions	first value	Recurrence steps
$z_{m}=f^{m}(c,z_{0})$	$z_{0}=0$	$z_{1}=c$	$z_{m+1}=z_{m}^{2}+c$
$P_{m}=f^{m}(c,z)-re^{2\pi i\theta }$			$P_{m}=z_{m}-t$
$D_{m}={\dfrac {\partial }{\partial c}}P_{m}={\dfrac {\partial }{\partial c}}f^{m}$	$D_{0}=0$	$D_{1}=1$	$D_{m+1}=2z_{m}D_{m}+1$

$t=re^{2\pi i\theta }=const$

Newton map :

arbitrary names	Initial conditions	first value	Recurrence steps
$N_{m+1}(c_{m})=c_{m}-{\frac {P_{m}}{D_{m}}}$			$N_{m+1}=c_{m}-{\frac {z_{m}-t}{D_{m}}}$

$c_{0}\quad {\xrightarrow {N}}\ c_{1}\quad {\xrightarrow {N}}\ c_{2}\quad {\xrightarrow {N}}\ c_{3}...$

Compute n points of the ray

 // compute and send to the output 4*depth complex points ( pair of arbitrary precision numbers )
  for (int i = 0; i < depth * 4; ++i) {
  
    mandelbrot_external_ray_in_step(ray);
    mandelbrot_external_ray_in_get(ray, cre, cim);
    // send new c to output  : (input of other program) or (text file  )
    mpfr_out_str(0, 10, 0, cre, GMP_RNDN); putchar(' ');
    mpfr_out_str(0, 10, 0, cim, GMP_RNDN); putchar('\n');
  }

Radial parameter r

Depth

Using fixed integer D ( depth ) : ^[21]

$d_{max}=D>1$

and fixed maximal value of radial parameter :

$r_{min}=R>1$

one can compute D points of ray using fomula :

$r=R^{1/2^{d}}$

which is :

$1<R^{1/2^{D}}<=r<=R$

When $d\to d_{max}=D\$ then $r\to 1$ and radius reaches enough close to the boundary of Mandelbrot set

/*
Maxima CAS code 

Number of ray points = depth  
r = radial parametr : 1 < R^{1/{2^D}} <= r > er  
*/


GiveRay( angle, depth):=
block (
 [ r, R: 65536],
 r:R,
 for d:1   thru depth  do
   (
     r:er^(1/(2^d)),
     print("d = ", d, " r = ", float(r))
    )
)$

compile(all)$


/* --------------------------- */

GiveRay(1/3,25)$

Output :

(%i4) GiveRay(1/3,25)
"d = "1" r = "256.0
"d = "2" r = "16.0
"d = "3" r = "4.0
"d = "4" r = "2.0
"d = "5" r = "1.414213562373095
"d = "6" r = "1.189207115002721
"d = "7" r = "1.090507732665258
"d = "8" r = "1.044273782427414
"d = "9" r = "1.021897148654117
"d = "10" r = "1.0108892860517
"d = "11" r = "1.005429901112803
"d = "12" r = "1.002711275050203
"d = "13" r = "1.001354719892108
"d = "14" r = "1.000677130693066
"d = "15" r = "1.000338508052682
"d = "16" r = "1.000169239705302
"d = "17" r = "1.000084616272694
"d = "18" r = "1.000042307241396
"d = "19" r = "1.000021153396965
"d = "20" r = "1.00001057664255
"d = "21" r = "1.000005288307292
"d = "22" r = "1.00000264415015
"d = "23" r = "1.000001322074201
"d = "24" r = "1.000000661036882
"d = "25" r = "1.000000330518386

Depth and sharpness

Using fixed integer S called sharpness :

$S>0$

and ...( to do )

one can compute S*D points of ray using fomula :

... ( to do)

Note that last point is the same as in depth method, but there are more points here ( S*D > D ) .

GiveRay( angle, depth, sharpness):=
block (
 [ r, R: 65536, m ],
 
 r:R,
 
 for d:1   thru depth  do
   (
     for s:1   thru sharpness  do
      (  m: (d-1)*sharpness +s,
         r:R^(1/(2^(m/sharpness))),
         print("d = ", d, " ; s = ", s , "; r = ", float(r))
      )
    )
)$

compile(all)$


/* --------------------------- */

GiveRay(1/3,25, 2)$