Fractals/Mathematics/Newton method

From Wikibooks, open books for an open world
Jump to: navigation, search

Introduction[edit]

types[edit]

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
  • a complex-valued function
  • a system of equations

What is needed ?[edit]

  • 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 )[3]


stopping rule[edit]

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


Pitfalls[edit]

  • method has a cycle and never converge
    • inflection point [4]
    • 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

Pseudocode[edit]

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[edit]

/* 
 
Maxima CAS code 
from /maxima/share/numeric/newton1.mac 
input : 
 exp  =  function of one variable, x
 var = variable 
 x0 = inital value of variable
 eps = 

The search begins with x = x_0 and proceeds until abs(expr) < eps (with expr evaluated at the current value of x). 

output : xn = an approximate solution of expr = 0 by Newton's method
*/ 


newton(exp,var,x0,eps):=
block(
  [xn,s,numer],
   numer:true,
   s:diff(exp,var),
   xn:x0,
   
   loop, 
    if abs(subst(xn,var,exp))<eps 
         then return(xn),
    xn:xn-subst(xn,var,exp)/subst(xn,var,s),
   go(loop) 
)$


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

recurrence relations[edit]

Define the iterated function f^n by:

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

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} = 2 A_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} = 2 A_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{ definition of f} \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{definition of A} 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} 2 A_m (\dfrac{\partial}{\partial z} A_m)  \quad \xrightarrow{definition\ of\ B} 2 A_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} (2 A_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} 2 A_m (\dfrac{\partial}{\partial c} A_m) + 1  \quad \xrightarrow{definition\ of\ D} 2 A_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} (2 A_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{definition of\ 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[edit]

A center n of period p satisfies :

f^p (0, n) = 0


Applying Newton's method in one variable:

\begin{align}
n_0 &= \text{ initial estimate } \\
n_{m+1} &= 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)}
\end{align}

boundary[edit]

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 i t} &= 0
\end{cases}

Defining a function of two complex variables:

\begin{align}
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 i t} \end{pmatrix}
\end{align}

and applying Newton's method in two variables:

\begin{align}
\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{align}

where

\begin{align}
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 i t}) &
\dfrac{\partial}{\partial c} (\dfrac{\partial }{\partial z}f^p(z, c) - e^{2 \pi i t})
\end{pmatrix}
\end{align}

This can be expressed using the recurrence relations as:

\begin{align}
\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 i t} \end{pmatrix}
\end{align}

where \{A,B,C,D,E\}_p are evaluated at (w_m, b_m).

Example[edit]

"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 [5])

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{align}
\boldsymbol{v}_0 = 
\begin{pmatrix} 
n\\
n
 \end{pmatrix}
=
\begin{pmatrix} 
-0.12256116687665+0.74486176661974i\\
-0.12256116687665+0.74486176661974i 
 \end{pmatrix}

\end{align}



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 i 0} &= 0
\end{cases}

so :


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

Then function g

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


Jacobian matrix J is


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


J^{-1}_g(\boldsymbol{v}) = 
\begin{pmatrix}
\frac{24z^5+48cz^3+8z^3+24c^2z+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^2z+16cz)(8z^7+24cz^5+24c^2z^3+8cz^3+8c^3z+8c^2z-1)+(-56z^6-120cz^4-72c^2z^2-24cz^2-8c^3-8c^2)(4z^6+12cz^4+2z^4+12c^2z^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^{-1}_g(\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^{-1}_g(\boldsymbol{v}_0)

 \boldsymbol{v}_2 = (\boldsymbol{v}_1 -g(\boldsymbol{v}_1) ) J^{-1}_g(\boldsymbol{v}_1)

 ...

 \boldsymbol{v}_{m+1} = (\boldsymbol{v}_m -g(\boldsymbol{v}_m) ) J^{-1}_g(\boldsymbol{v}_m)

Maxima CAS code[edit]

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

In Maxima CAS one can use Newton method for solving systems of multiple nonlinear functions. It is implemented in mnewton function. See directory :

/Maxima..../share/contrib/mnewton.mac

which uses linsolve_by_lu defined in :

 /share/linearalgebra/lu.lisp

See this image for more code.

C++ code[edit]

/*
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[edit]

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:

  1. given a center location estimate n_m accurate to P bits of precision
  2. compute a boundary location estimate b_m accurate to the same precision, using center n_m as starting value
  3. 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.

external rays[edit]

Dynamic external ray[edit]

//  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[edit]

External parameter rays

Compute one point of the ray[edit]

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) - r e^{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} = 2 z_m D_m + 1


t =  r e^{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[edit]

 // 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[edit]

Depth


Using fixed integer D ( depth ) : [6]

 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)$

Output :

(%i4) GiveRay(1/3,25,2)
"d = "1" ; s = "1"; r = "2545.456152628088
"d = "1" ; s = "2"; r = "256.0
"d = "2" ; s = "1"; r = "50.45251383854013
"d = "2" ; s = "2"; r = "16.0
"d = "3" ; s = "1"; r = "7.10299330131601
"d = "3" ; s = "2"; r = "4.0
"d = "4" ; s = "1"; r = "2.665144142690224
"d = "4" ; s = "2"; r = "2.0
"d = "5" ; s = "1"; r = "1.632526919438152
"d = "5" ; s = "2"; r = "1.414213562373095
"d = "6" ; s = "1"; r = "1.277703768264832
"d = "6" ; s = "2"; r = "1.189207115002721
"d = "7" ; s = "1"; r = "1.13035559372475
"d = "7" ; s = "2"; r = "1.090507732665258
"d = "8" ; s = "1"; r = "1.063181825335982
"d = "8" ; s = "2"; r = "1.044273782427414
"d = "9" ; s = "1"; r = "1.031107087230023
"d = "9" ; s = "2"; r = "1.021897148654117
"d = "10" ; s = "1"; r = "1.015434432757735
"d = "10" ; s = "2"; r = "1.0108892860517
"d = "11" ; s = "1"; r = "1.007687666272509
"d = "11" ; s = "2"; r = "1.005429901112803
"d = "12" ; s = "1"; r = "1.003836473870375
"d = "12" ; s = "2"; r = "1.002711275050203
"d = "13" ; s = "1"; r = "1.001916400639482
"d = "13" ; s = "2"; r = "1.001354719892108
"d = "14" ; s = "1"; r = "1.000957741685173
"d = "14" ; s = "2"; r = "1.000677130693066
"d = "15" ; s = "1"; r = "1.000478756238819
"d = "15" ; s = "2"; r = "1.000338508052682
"d = "16" ; s = "1"; r = "1.000239349475324
"d = "16" ; s = "2"; r = "1.000169239705302
"d = "17" ; s = "1"; r = "1.000119667577497
"d = "17" ; s = "2"; r = "1.000084616272694
"d = "18" ; s = "1"; r = "1.000059831998815
"d = "18" ; s = "2"; r = "1.000042307241396
"d = "19" ; s = "1"; r = "1.000029915551937
"d = "19" ; s = "2"; r = "1.000021153396965
"d = "20" ; s = "1"; r = "1.000014957664103
"d = "20" ; s = "2"; r = "1.00001057664255
"d = "21" ; s = "1"; r = "1.000007478804085
"d = "21" ; s = "2"; r = "1.000005288307292
"d = "22" ; s = "1"; r = "1.000003739395051
"d = "22" ; s = "2"; r = "1.00000264415015
"d = "23" ; s = "1"; r = "1.000001869695778
"d = "23" ; s = "2"; r = "1.000001322074201
"d = "24" ; s = "1"; r = "1.000000934847452
"d = "24" ; s = "2"; r = "1.000000661036882
"d = "25" ; s = "1"; r = "1.000000467423617
"d = "25" ; s = "2"; r = "1.000000330518386

Src code[edit]

/* 
 
it is a c console program which 
- computes external parameter ray 
- for complex quadratic polynomial fc(z) = z^2+c
- uses arbitrary precision ( mpfr) with dynamic precision adjustment
- uses Newton method ( described by T Kawahira) http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf
 
 
 
gcc  -o rayin e.c -std=c99 -Wall -Wextra -pedantic -lmpfr -lgmp -lm 
 
 ./rayin 1/3 25
 ./rayin 1/3 25 >ray13.txt
 ./rayin 1/3 25  | ./rescale 53 53 -0.75 0 1.5 0 >rray13.txt
 
 
 
 program uses the code by Claude Heiland-Allen
 from https://gitorious.org/maximus/book/
 
 
 
 
" ray_in computes 4 points per dwell band, and ray_out currently computes 16.  
 Moreover ray_in has dynamic precision adjustment and 
ray_out doesn't (so you need a high precision to get through the thin 
gaps) - fixing both of these is on my mental todo list.  ray_out will 
always be slower though, as it needs to compute extra points when 
crossing dwell bands to accumulate angle bits (and to work correctly)."
 
What algorithm do you use for drawing external ray ?
 
"essentially the Newton's method one in
http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf
 
precision is automatic, effectively it checks how close dwell bands are 
together and uses higher precision when they are narrow and lower 
precision when they are wide.  in general precision will increase along 
the ray as it gets closer to the boundary. "
 
 
"periodic rays stop near the landing point at the current zoom level, 
manually retrace if you zoom in and find the gap annoying.  tracing 
periodic rays is a lot slower than tracing preperiodic rays, especially 
if you zoom in to the cusp - suggestions welcome for improvements here 
(perhaps a larger radius than a single pixel for landing point nearness 
threshold?). "
 
 
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h> // arbitrary precision 
 
 
 
 
const float pif = 3.14159265358979323846264338327950288419716939937510f;
const double pi = 3.14159265358979323846264338327950288419716939937510;
const long double pil = 3.14159265358979323846264338327950288419716939937510l;
 
 
 
struct mandelbrot_external_ray_in {
 
  mpq_t angle;
  mpq_t one;
  unsigned int sharpness; // number of steps to take within each dwell band
  unsigned int precision; // delta needs this many bits of effective precision
  unsigned int accuracy;  // epsilon is scaled relative to magnitude of delta
  double escape_radius;
  mpfr_t epsilon2;
  mpfr_t cx;
  mpfr_t cy;
  unsigned int j;
  unsigned int k;
  // temporaries
  mpfr_t zx, zy, dzx, dzy, ux, uy, vx, vy, ccx, ccy, cccx, cccy;
};
 
 
 
// new : allocates memory, inits and sets variables
// input angle = external angle in turns 
extern struct mandelbrot_external_ray_in *mandelbrot_external_ray_in_new(mpq_t angle) {
 
  struct mandelbrot_external_ray_in *r = malloc(sizeof(struct mandelbrot_external_ray_in));
 
  mpq_init(r->angle);
  mpq_set(r->angle, angle);
 
  mpq_init(r->one);
  mpq_set_ui(r->one, 1, 1);
 
  r->sharpness = 4;
  r->precision = 4;
  r->accuracy  = 4;
  r->escape_radius = 65536.0;
 
  mpfr_init2(r->epsilon2, 53);
  mpfr_set_ui(r->epsilon2, 1, GMP_RNDN);
  // external angle in radions
  double a = 2.0 * pi * mpq_get_d(r->angle);
  // initial value c = 
  mpfr_init2(r->cx, 53);
  mpfr_init2(r->cy, 53);
  mpfr_set_d(r->cx, r->escape_radius * cos(a), GMP_RNDN);
  mpfr_set_d(r->cy, r->escape_radius * sin(a), GMP_RNDN);
 
  r->k = 0;
  r->j = 0;
  // initialize temporaries
  mpfr_inits2(53, r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0);
  return r;
}
 
 
// delete
extern void mandelbrot_external_ray_in_delete(struct mandelbrot_external_ray_in *r) {
  mpfr_clear(r->epsilon2);
  mpfr_clear(r->cx);
  mpfr_clear(r->cy);
  mpq_clear(r->angle);
  mpq_clear(r->one);
  mpfr_clears(r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0);
  free(r);
}
 
 
 
// step 
extern int mandelbrot_external_ray_in_step(struct mandelbrot_external_ray_in *r) {
 
  if (r->j >= r->sharpness) {
    mpq_mul_2exp(r->angle, r->angle, 1);
    if (mpq_cmp_ui(r->angle, 1, 1) >= 0) {
      mpq_sub(r->angle, r->angle, r->one);
    }
    r->k++;
    r->j = 0;
  }
  // r0 <- er ** ((1/2) ** ((j + 0.5)/sharpness))
  // t0 <- r0 cis(2 * pi * angle)
  double r0 = pow(r->escape_radius, pow(0.5, (r->j + 0.5) / r->sharpness));
  double a0 = 2.0 * pi * mpq_get_d(r->angle);
  double t0x = r0 * cos(a0);
  double t0y = r0 * sin(a0);
  // c <- r->c
  mpfr_set(r->ccx, r->cx, GMP_RNDN);
  mpfr_set(r->ccy, r->cy, GMP_RNDN);
  for (unsigned int i = 0; i < 64; ++i) { // FIXME arbitrary limit
    // z <- 0
    // dz <- 0
    mpfr_set_ui(r->zx, 0, GMP_RNDN);
    mpfr_set_ui(r->zy, 0, GMP_RNDN);
    mpfr_set_ui(r->dzx, 0, GMP_RNDN);
    mpfr_set_ui(r->dzy, 0, GMP_RNDN);
    // iterate
    for (unsigned int p = 0; p <= r->k; ++p) {
      // dz <- 2 z dz + 1
      mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN);
      mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN);
      mpfr_mul(r->vx, r->zx, r->dzy, GMP_RNDN);
      mpfr_mul(r->vy, r->zy, r->dzx, GMP_RNDN);
      mpfr_sub(r->dzx, r->ux, r->uy, GMP_RNDN);
      mpfr_add(r->dzy, r->vx, r->vy, GMP_RNDN);
      mpfr_mul_2ui(r->dzx, r->dzx, 1, GMP_RNDN);
      mpfr_mul_2ui(r->dzy, r->dzy, 1, GMP_RNDN);
      mpfr_add_ui(r->dzx, r->dzx, 1, GMP_RNDN);
      // z <- z^2 + c
      mpfr_sqr(r->ux, r->zx, GMP_RNDN);
      mpfr_sqr(r->uy, r->zy, GMP_RNDN);
      mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN);
      mpfr_mul(r->vy, r->zx, r->zy, GMP_RNDN);
      mpfr_mul_2ui(r->vy, r->vy, 1, GMP_RNDN);
      mpfr_add(r->zx, r->vx, r->ccx, GMP_RNDN);
      mpfr_add(r->zy, r->vy, r->ccy, GMP_RNDN);
    }
    // c' <- c - (z - t0) / dz
    mpfr_sqr(r->ux, r->dzx, GMP_RNDN);
    mpfr_sqr(r->uy, r->dzy, GMP_RNDN);
    mpfr_add(r->vy, r->ux, r->uy, GMP_RNDN);
    mpfr_sub_d(r->zx, r->zx, t0x, GMP_RNDN);
    mpfr_sub_d(r->zy, r->zy, t0y, GMP_RNDN);
    mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN);
    mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN);
    mpfr_add(r->vx, r->ux, r->uy, GMP_RNDN);
    mpfr_div(r->ux, r->vx, r->vy, GMP_RNDN);
    mpfr_sub(r->cccx, r->ccx, r->ux, GMP_RNDN);
    mpfr_mul(r->ux, r->zy, r->dzx, GMP_RNDN);
    mpfr_mul(r->uy, r->zx, r->dzy, GMP_RNDN);
    mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN);
    mpfr_div(r->uy, r->vx, r->vy, GMP_RNDN);
    mpfr_sub(r->cccy, r->ccy, r->uy, GMP_RNDN);
    // delta^2 = |c' - c|^2
    mpfr_sub(r->ux, r->cccx, r->ccx, GMP_RNDN);
    mpfr_sub(r->uy, r->cccy, r->ccy, GMP_RNDN);
    mpfr_sqr(r->vx, r->ux, GMP_RNDN);
    mpfr_sqr(r->vy, r->uy, GMP_RNDN);
    mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN);
    // enough_bits = 0 < 2 * prec(eps^2) + exponent eps^2 - 2 * precision
    int enough_bits = 0 < 2 * (mpfr_get_prec(r->epsilon2) - r->precision) + mpfr_get_exp(r->epsilon2);
    if (enough_bits) {
      // converged = delta^2 < eps^2
      int converged = mpfr_less_p(r->ux, r->epsilon2);
      if (converged) {
        // eps^2 <- |c' - r->c|^2 >> (2 * accuracy)
        mpfr_sub(r->ux, r->cccx, r->cx, GMP_RNDN);
        mpfr_sub(r->uy, r->cccy, r->cy, GMP_RNDN);
        mpfr_sqr(r->vx, r->ux, GMP_RNDN);
        mpfr_sqr(r->vy, r->uy, GMP_RNDN);
        mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN);
        mpfr_div_2ui(r->epsilon2, r->ux, 2 * r->accuracy, GMP_RNDN);
        // j <- j + 1
        r->j = r->j + 1;
        // r->c <- c'
        mpfr_set(r->cx, r->cccx, GMP_RNDN);
        mpfr_set(r->cy, r->cccy, GMP_RNDN);
        return 1;
      }
    } else {
      // bump precision
      mpfr_prec_t prec = mpfr_get_prec(r->cx) + 32;
      mpfr_prec_round(r->cx, prec, GMP_RNDN);
      mpfr_prec_round(r->cy, prec, GMP_RNDN);
      mpfr_prec_round(r->epsilon2, prec, GMP_RNDN);
      mpfr_set_prec(r->ccx, prec);
      mpfr_set_prec(r->ccy, prec);
      mpfr_set_prec(r->cccx, prec);
      mpfr_set_prec(r->cccy, prec);
      mpfr_set_prec(r->zx, prec);
      mpfr_set_prec(r->zy, prec);
      mpfr_set_prec(r->dzx, prec);
      mpfr_set_prec(r->dzy, prec);
      mpfr_set_prec(r->ux, prec);
      mpfr_set_prec(r->uy, prec);
      mpfr_set_prec(r->vx, prec);
      mpfr_set_prec(r->vy, prec);
      i = 0;
    }
    mpfr_set(r->ccx, r->cccx, GMP_RNDN);
    mpfr_set(r->ccy, r->cccy, GMP_RNDN);
  }
  return 0;
}
 
 
 
// get 
extern void mandelbrot_external_ray_in_get(struct mandelbrot_external_ray_in *r, mpfr_t x, mpfr_t y) {
  mpfr_set_prec(x, mpfr_get_prec(r->cx));
  mpfr_set(x, r->cx, GMP_RNDN);
  mpfr_set_prec(y, mpfr_get_prec(r->cy));
  mpfr_set(y, r->cy, GMP_RNDN);
}
 
 
 
 
void usage(const char *progname) {
  fprintf(stderr,
    "usage: %s angle depth\n"
    "outputs space-separated complex numbers on stdout\n"
    "example : \n %s 1/3 25\n"
    "or \n %s 1/3 25 > r.txt\n",
    progname, progname, progname);
}
 
 
 
 
 
// needs input 
int main(int argc, char **argv) {
 
  // read and check input
  if (argc < 3) { usage(argv[0]); return 1; }
 
 
  mpq_t angle;
  mpq_init(angle);
  mpq_set_str(angle, argv[1], 0);
  mpq_canonicalize(angle);
 
  int depth = atoi(argv[2]);
 
  struct mandelbrot_external_ray_in *ray = mandelbrot_external_ray_in_new(angle);
 
  mpfr_t cre, cim;
  mpfr_init2(cre, 53);
  mpfr_init2(cim, 53);
 
  // compute and send output
  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  
    mpfr_out_str(0, 10, 0, cre, GMP_RNDN); putchar(' ');
    mpfr_out_str(0, 10, 0, cim, GMP_RNDN); putchar('\n');
  }
 
 
 
  // clear 
  mpfr_clear(cre);
  mpfr_clear(cim);
  mandelbrot_external_ray_in_delete(ray);
  mpq_clear(angle);
 
 return 0;
}

See also[edit]

References[edit]

  1. NEWTON’S METHOD AND FRACTALS by AARON BURTON
  2. Newton, Chebyshev, and Halley Basins of Attraction; A Complete Geometric Approach by Bart D. Stewart
  3. Stopping Criterion by Prof. Alan Hood
  4. Numerical Methods for Scientists and Engineers by Richard Hamming, page 69
  5. Home page of Claude Heiland-Allen
  6. Drawing external ray using Newton method ( described by T Kawahira)