Fractals/Mathematics/Newton method

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

Introduction

recurrence relations

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 derivatives:

$\begin{align} A_m &= f^m \\ B_m &= \dfrac{\partial}{\partial z} f^m \\ C_m &= \dfrac{\partial}{\partial z} \dfrac{\partial}{\partial z} f^m \\ D_m &= \dfrac{\partial}{\partial c} f^m \\ E_m &= \dfrac{\partial}{\partial c} \dfrac{\partial}{\partial z} f^m \end{align}$

These derivatives can be computed by recurrence relations.

Initial conditions:

$\begin{align} A_0 &= z \\ B_0 &= 1 \\ C_0 &= 0 \\ D_0 &= 0 \\ E_0 &= 0 \end{align}$

Recurrence steps:

$\begin{align} A_{m+1} &= A_m ^ 2 + c \\ B_{m+1} &= 2 A_m B_m \\ C_{m+1} &= 2 (B_m ^ 2 + A_m C_m) \\ D_{m+1} &= 2 A_m D_m + 1 \\ E_{m+1} &= 2 (A_m E_m + B_m D_m) \end{align}$

derivation of recurrence relations

Derivation of $A_0$:

$\begin{align} & A_0 & \\ = & & \text{ ( definition of }A\text{ ) } \\ & f^0 (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & z & \end{align}$

Derivation of $B_0$:

$\begin{align} & B_0 & \\ = & & \text{ ( definition of }B\text{ ) } \\ & \dfrac{\partial}{\partial z} f^0 (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & \dfrac{\partial}{\partial z} z & \\ = & & \text{ ( derivative ) } \\ & 1 & \end{align}$

Derivation of $C_0$:

$\begin{align} & C_0 & \\ = & & \text{ ( definition of }C\text{ ) } \\ & \dfrac{\partial}{\partial z} \dfrac{\partial}{\partial z} f^0 (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & \dfrac{\partial}{\partial z} \dfrac{\partial}{\partial z} z & \\ = & & \text{ ( derivative ) } \\ & \dfrac{\partial}{\partial z} 1 & \\ = & & \text{ ( derivative ) } \\ & 0 & \end{align}$

Derivation of $D_0$:

$\begin{align} & D_0 & \\ = & & \text{ ( definition of }D\text{ ) } \\ & \dfrac{\partial}{\partial c} f^0 (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & \dfrac{\partial}{\partial c} z & \\ = & & \text{ ( derivative ) } \\ & 0 & \end{align}$

Derivation of $E_0$:

$\begin{align} & E_0 & \\ = & & \text{ ( definition of }E\text{ ) } \\ & \dfrac{\partial}{\partial c} \dfrac{\partial}{\partial z} f^0 (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & \dfrac{\partial}{\partial c} \dfrac{\partial}{\partial z} z & \\ = & & \text{ ( derivative ) } \\ & \dfrac{\partial}{\partial c} 1 & \\ = & & \text{ ( derivative ) } \\ & 0 & \end{align}$

Derivation of $A_{m+1}$:

$\begin{align} & A_{m+1} & \\ = & & \text{ ( definition of }A\text{ ) } \\ & f^{m+1} (z, c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & f (f^m (z, c), c) & \\ = & & \text{ ( definition of }A\text{ ) } \\ & f (A_m , c) & \\ = & & \text{ ( definition of }f\text{ ) } \\ & A_m ^ 2 + c & \end{align}$

Derivation of $B_{m+1}$:

$\begin{align} & B_{m+1} & \\ = & & \text{ ( definition of }B\text{ ) } \\ & \dfrac{\partial}{\partial z} A_{m+1} & \\ = & & \text{ ( definition of }A\text{ ) } \\ & \dfrac{\partial}{\partial z} (A_m ^ 2 + c) & \\ = & & \text{ ( distributivity ) } \\ & \dfrac{\partial}{\partial z} A_m ^ 2 + \dfrac{\partial}{\partial z} c & \\ = & & \text{ ( constant derivative ) } \\ & \dfrac{\partial}{\partial z} A_m ^ 2 + 0 & \\ = & & \text{ ( zero ) } \\ & \dfrac{\partial}{\partial z} A_m ^ 2 & \\ = & & \text{ ( chain rule ) } \\ & 2 A_m (\dfrac{\partial}{\partial z} A_m) & \\ = & & \text{ ( definition of }B\text{ ) } \\ & 2 A_m B_m & \end{align}$

Derivation of $C_{m+1}$:

$\begin{align} & C_{m+1} & \\ = & & \text{ ( definition of }C\text{ ) } \\ & \dfrac{\partial}{\partial z} B_{m+1} & \\ = & & \text{ ( definition of }B\text{ ) } \\ & \dfrac{\partial}{\partial z} (2 A_m B_m) & \\ = & & \text{ ( linearity ) } \\ & 2 \dfrac{\partial}{\partial z} (A_m B_m) & \\ = & & \text{ ( product rule ) } \\ & 2 ( (\dfrac{\partial}{\partial z} A_m) B_m + A_m (\dfrac{\partial}{\partial z} B_m) ) & \\ = & & \text{ ( definition of }B\text{ ) } \\ & 2 ( B_m B_m + A_m (\dfrac{\partial}{\partial z} B_m) ) & \\ = & & \text{ ( algebra ) } \\ & 2 ( B_m ^ 2 + A_m (\dfrac{\partial}{\partial z} B_m) ) & \\ = & & \text{ ( definition of }C\text{ ) } \\ & 2 ( B_m ^ 2 + A_m C_m ) & \end{align}$

Derivation of $D_{m+1}$:

$\begin{align} & D_{m+1} & \\ = & & \text{ ( definition of }D\text{ ) } \\ & \dfrac{\partial}{\partial c} A_{m+1} & \\ = & & \text{ ( definition of }A\text{ ) } \\ & \dfrac{\partial}{\partial c} (A_m ^ 2 + c) & \\ = & & \text{ ( distributivity ) } \\ & \dfrac{\partial}{\partial c} A_m ^ 2 + \dfrac{\partial}{\partial c} c & \\ = & & \text{ ( derivative ) } \\ & \dfrac{\partial}{\partial c} A_m ^ 2 + 1 & \\ = & & \text{ ( chain rule ) } \\ & 2 A_m (\dfrac{\partial}{\partial c} A_m) + 1 & \\ = & & \text{ ( definition of }D\text{ ) } \\ & 2 A_m D_m + 1 & \end{align}$

Derivation of $E_{m+1}$:

$\begin{align} & E_{m+1} & \\ = & & \text{ ( definition of }E\text{ ) } \\ & \dfrac{\partial}{\partial c} B_{m+1} & \\ = & & \text{ ( definition of }B\text{ ) } \\ & \dfrac{\partial}{\partial c} (2 A_m B_m) & \\ = & & \text{ ( linearity ) } \\ & 2 \dfrac{\partial}{\partial c} (A_m B_m) & \\ = & & \text{ ( product rule ) } \\ & 2 ( A_m (\dfrac{\partial}{\partial c} B_m) + (\dfrac{\partial}{\partial c} A_m) B_m ) & \\ = & & \text{ ( definition of }E\text{ ) } \\ & 2 ( A_m E_m + (\dfrac{\partial}{\partial c} A_m) B_m ) & \\ = & & \text{ ( definition of }D\text{ ) } \\ & 2 ( A_m E_m + D_m B_m ) & \end{align}$

center

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)} \end{align}$

or equivalently:

$\begin{align} n_0 &= \text{ initial estimate } \\ n_{m+1} &= n_m - \frac{A_p (0, n_m)}{D_p (0, n_m)} \end{align}$

boundary

Boundaries of Components of Mandelbrot set by Newton method

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

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

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

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

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:

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

External rays

See also

• The Newton-Raphson fractal

References