Fractals/Mathematics/Derivative

• Derivative of Iterated function (map)
• of the function f with respect to variable ...
• following the derivative
• the rolling derivative
• the derivative of order n: first, second, ...
• Keep track of derivatives ...

notation

• variable z
• constant parameter c
• function f
• derivative d
• iteration number p
• wrt = with respect to

Notation types

• Leibniz's notation
• Euler's notation
• Newton's notation ( dot notation)
• Lagrange's notation ( prime mark notation)

${\frac {d}{dz}}f(z)=f'(z)=z'=d=D\,$

Informal description

The derivative gives you what is called the gradient. For 2D, consider the terrain as a contour map. The gradient tells you the rate of change in x and y directions. The slope is the magnitude of the gradient, and the contours are closer together for higher slope. For 3D the gradient is the rate of change in x, y and z, and the slope is still the magnitude of the gradient. Technically, slope is only 1D, and the gradient also tells you the direction of the slope. When we talk about derivative giving you the slope we are being casual. ( xenodreambuie )

Derivative of the iterated function

  "the derivative basically as it's calculated for anlytical DE using the implementation of the chain rule for the derivative of f(g(x)) where f(x) is say g(x)^p+c and the value of g(x) is current z. So f'(g(x)) is p*z^(p-1)  and g'(x) is the derivative from the previous iteration,  so on each iteration the new derivative is:  dz = dz*p*z^(p-1)  and new z = z^p+c as normal" FractalDave

For Julia set:

// initial values:
c is const 
z =  z0 ( z is a variable, read initial value from the pixel)
dz = 1


// Recurrence steps: 
dz = f'(z)*dz   // calculating dz and multiply the previous d
z = f(z)       // forward iteration

For Mandelbrot set:

// initial values:
z = critical point ( critical pointis  0 for multibrot sets, z is variable)
c = ( c is a variable, read initial value from the pixel)
dc = 0


// Recurrence steps: 
dc = f'(z)*dc  + 1 // calculating dc and multiply the previous dc
z = f(z)       // forward iteration

It can be computed with Maxima CAS :

display2d:false;
f:1/(z^3+d*z);
dz : diff(f,z,1); // first derivative of f wrt variable z 
dc :  diff(f,c,1); // first derivative of f wrt parameter c

Rational function

 $f(x)={\frac {P(x)}{Q(x)}}$

 $f'(x)={\frac {P'(x)Q(x)-P(x)Q'(x)}{Q^{2}(x)}}$

degree 3

wrt z

// function 
f(z)= 1/(z^3 + c*z )

// first derivativwe wrt z
d = f'(z) = -(3z^2 +c)/(z^3 + cz)^2

// iteration:
z_(n+1) = f(z_n)


// initial values:
z = z
d = 1

// Recurrence steps:
d = - d*(3z^2 +c)/(z^3 + cz)^2
z = 1/(z^3 + c*z) 

degree 5

Julia set f(z)=1 over az5+z3+bz

display2d:false;
f: 1/(a*z^5+ z^3 + b*z);
dz: diff(f,z,1);

in plain text:

  dz = -(5*a*z^4+3*z^2+b)/(a*z^5+z^3+b*z)^2

Example:

• 1/((0.15*I+0.15)*z^5+(3*I-3)*z + z^3)

degree 5 by Buschmann

Julia set drawn by distance estimation, the iteration is of the form $f_{c}(z)=1-z^{2}+{\frac {z^{5}}{2+4z}}+c$ in black and white

f(z):=z^5/(4*z+2)-z^2+c+1

(%i4) diff(f(z), z,1);

(%o4) (5*z^4)/(4*z+2)-(4*z^5)/(4*z+2)^2-2*z

Polynomials

First derivative wrt z

degree	Function f(z)	derivative wrt z
2	$z^{2}+c$	$2*z$
3	$z^{3}+c$	$3*z^{2}$
4	$z^{4}+c$	$4*z^{3}$
d	$z^{d}+c$	$d*z^{d-1}$

C src code for making exterior DEM/M images of Multibrot sets: Mandelbrot sets for $z^{q}+c$ where q is a degree of unicritical monomial:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> //OpenM

/*
fork of 
mandelbrot-book how to write a book about the Mandelbrot set by Claude Heiland-Alle
https://code.mathr.co.uk/mandelbrot-book/blob/HEAD:/book/

https://stackoverflow.com/questions/77135883/why-do-my-functions-not-work-in-parallel


gcc s.c -lm -Wall -fopenmp

./a.out >ed.ppm   // P6 = binary Portable PixMap see en wikipedia: Netpbm#File_formats


*/

const double pi = 3.141592653589793;

/*
 int q = 5 ;
complex double f(complex double z, int q){ return cpow(z,q) + c;}
complex double d(complex double z, int q) {return q*cpow(z, q-1); }

*/
complex double f(complex double z, complex double c){ return z*z*z*z*z + c;}
complex double d(complex double z) {return 5*z*z*z*z; }




double cnorm(double _Complex z) // https://stackoverflow.com/questions/6363247/what-is-a-complex-data-type-and-an-imaginary-data-type-in-c
{
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

void hsv2rgb(double h, double s, double v, int *red, int *grn, int *blu) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *red = fmin(fmax(255 * r + 0.5, 0), 255);
  *grn = fmin(fmax(255 * g + 0.5, 0), 255);
  *blu = fmin(fmax(255 * b + 0.5, 0), 255);
}

int main()
{
  const int aa = 4;
  const int w = 800 * aa;
  const int h = 800 * aa;
  const int n = 1024;
  const double r = 2;
  const double px = r / (h/2);
  const double r2 = 25 * 25;
  unsigned char *img = malloc(3 * w * h);

  #pragma omp parallel for schedule(dynamic)
  for (int j = 0; j < h; ++j)
  {
    double _Complex c;
    double _Complex z;
    double _Complex dc;
    double y = (h/2 - (j + 0.5)) / (h/2) * r;
    for (int i = 0; i < w; ++i)
    {
      double x =  (i + 0.5 - w/2) / (h/2) * r; // for q=2 add -0.5
      c = x + I * y;
      //double _Complex 
      dc = 0; // first derivative of zn with respect to c
      //double _Complex z = 0;
      z = 0;
      int k;
      for (k = 0; k < n; ++k)
      { 
      // 
        dc = d(z)*dc +1;
        z = f(z,c);
        
        if (cnorm(z) > r2)
          break;
      }
      
      // color
      double hue = 0, sat = 0, val = 1; // interior color = white
      
      if (k < n) 
      { // exterior and boundary color
        double _Complex de = 2 * z * log(cabs(z)) / dc;
        hue = fmod(1 + carg(de) / (2 * pi), 1); // ? slope of de
        sat = 0.25;
        val = tanh(cabs(de) / px / aa);
      }
      
      // hsv to rgb conversion
      int red, grn, blu;
      hsv2rgb(hue, sat, val, &red, &grn, &blu);
      // save rgb color to array
      img[3*(j * w + i)+0] = red;
      img[3*(j * w + i)+1] = grn;
      img[3*(j * w + i)+2] = blu;
    }
  }
  
  //
  printf("P6\n%d %d\n255\n", w, h);
  fwrite(img, 3 * w * h, 1, stdout);
  free(img);
  
  
  return 0;
}

C src code for making exterior DEM/J images of Julia sets for $z^{q}+c$ where q is a degree of unicritical monomial:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> //OpenM

/*
fork of 
mandelbrot-book how to write a book about the Mandelbrot set by Claude Heiland-Alle
https://code.mathr.co.uk/mandelbrot-book/blob/HEAD:/book/

https://stackoverflow.com/questions/77135883/why-do-my-functions-not-work-in-parallel


gcc j.c -lm -Wall -fopenmp

./a.out >ed.ppm   // P6 = binary Portable PixMap see en.wikipedia:  Netpbm#File_formats


*/



complex double examples [6] = {0, 0, I, -0.0649150006787816892861875745218343125883 +0.76821968591243610206311097043854440463 *I, -0.239026451915233+0.415516996006582*I, 0};

const double pi = 3.141592653589793;
int q = 4 ; // degree of multibrot set

/*
 
complex double f(complex double z, int q){ return cpow(z,q) + c;}
complex double d(complex double z, int q) {return q*cpow(z, q-1); }

*/
complex double f(complex double z, complex double c){ return z*z*z*z + c;} // multibrot z^q + c
complex double d(complex double z) {return 4*z*z*z; } // q*z^{q-1}




double cnorm(double _Complex z) // https://stackoverflow.com/questions/6363247/what-is-a-complex-data-type-and-an-imaginary-data-type-in-c
{
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

void hsv2rgb(double h, double s, double v, int *red, int *grn, int *blu) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *red = fmin(fmax(255 * r + 0.5, 0), 255);
  *grn = fmin(fmax(255 * g + 0.5, 0), 255);
  *blu = fmin(fmax(255 * b + 0.5, 0), 255);
}

int main()
{
  const int aa = 4;
  const int w = 800 * aa;
  const int h = 800 * aa;
  const int n = 1024;
  const double r = 2;
  const double px = r / (h/2);
  const double r2 = 25 * 25;
  unsigned char *img = malloc(3 * w * h);
  
  double _Complex c = examples[q];

  #pragma omp parallel for schedule(dynamic)
  for (int j = 0; j < h; ++j)
  {
    
        
    double y = (h/2 - (j + 0.5)) / (h/2) * r;
    for (int i = 0; i < w; ++i)
    {
      double x =  (i + 0.5 - w/2) / (h/2) * r; // for q=2 add -0.5
      double _Complex z = x + I * y;
      double _Complex dz = 1.0; // first derivative of zn with respect to z
      
      
      int k;
      for (k = 0; k < n; ++k)
      { 
        // 
        dz = d(z)*dz ;
        z = f(z,c);
        
        if (cnorm(z) > r2)
          break;
      }
      
      // color
      double hue = 0, sat = 0, val = 1; // interior color = white
      
      if (k < n) 
      { // exterior and boundary color
        double _Complex de = 2 * z * log(cabs(z)) / dz;
        hue = fmod(1 + carg(de) / (2 * pi), 1); // ? slope of de
        sat = 0.25;
        val = tanh(cabs(de) / px / aa);
      }
      
      // hsv to rgb conversion
      int red, grn, blu;
      hsv2rgb(hue, sat, val, &red, &grn, &blu);
      // save rgb color to array
      img[3*(j * w + i)+0] = red;
      img[3*(j * w + i)+1] = grn;
      img[3*(j * w + i)+2] = blu;
    }
  }
  
  //
  printf("P6\n%d %d\n255\n", w, h);
  fwrite(img, 3 * w * h, 1, stdout);
  free(img);
  
  
  return 0;
}

• 

Complex quadratic polynomial

First derivative wrt c

On parameter plane :

• $c$ is a variable
• $z_{0}=0$ is constant

${\frac {d}{dc}}f_{c}^{(p)}(z_{0})=z'_{p}\,$

This derivative can be found by iteration starting with

$z_{0}=0\,$

$z'_{0}=1\,$

and then ( compute derivative before next z because it uses previous values of z):

$z'_{p+1}=2\cdot z_{p}\cdot z'_{p}+1\,$

$z_{p+1}=z_{p}^{2}+c\,$

This can be verified by using the chain rule for the derivative.

${\frac {dz_{n+1}}{dc}}={\frac {d(z_{n}^{2}+c)}{dc}}=2z_{n}{\frac {dz_{n}}{dc}}+{\frac {dc}{dc}}=2z_{n}{\frac {dz_{n}}{dc}}+1\,$

• Maxima CAS function :

dcfn(p, z, c) :=
  if p=0 then 1
  else 2*fn(p-1,z,c)*dcfn(p-1, z, c)+1;

Example values :

$z_{0}=0\qquad \qquad z'_{0}=1\,$

$z_{1}=c\qquad \qquad z'_{1}=1\,$

$z_{2}=c^{2}+c\qquad z'_{2}=2c+1\,$

It can be used for:

• the distance estimation method for drawing a Mandelbrot set ( DEM/M )

• Mandelbrot set boundary = DEM/M
• 

First derivative wrt z

$z'_{n}\,$ is first derivative with respect to z.

This derivative can be found by iteration starting with

$z_{0}=z\,$

$z'_{0}=1\,$

and then :

$z'_{n+1}=2*z_{n}*z'_{n}\,$

$z_{n+1}=z_{n}^{2}+c\,$

description	arbitrary name	formula	Initial conditions	Recurrence step	common names
iterated complex quadratic polynomial	$A_{m}$	$f^{m}$	$A_{0}=z$	$A_{m+1}=A_{m}^{2}+c$	z or f
first derivative with respect to z	$B_{m}$	${\dfrac {\partial }{\partial z}}f^{m}$	$B_{0}=1$	$B_{m+1}=2A_{m}B_{m}$	dz, d (from derivative) or p ( from prime) of f'

Derivation of recurrence relation:

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

$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_{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$

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

Following the derivative

"As we iterate z, we can look at the derivatives of P at z. In our case it is quite simple: P′(z)=2z. Multiplying all these numbers along an orbit yields the derivative at z of the composition Pn. This multiplication can be carried on iteratively as we iterate z " A Cheritat

 der = 1
 z = c # note that we start with c instead of 0, to avoid multiplying the derivative by 0
 for n in range(0,N):
   new_z = z*z+c
   new_der = der*2*z
   z = new_z
   der = new_der

It can be used for :

• computing the external distance to the Julia set (DEM/J)
• detection of interior points of Mandelbrot set componnets

unsigned char ComputeColorOfDEMJ(complex double z){
	


	int nMax = IterMax_DEM;
	complex double dz = 1.0; //  is first derivative with respect to z.
	double distance;
	double cabsz;
	
	int n;

	for (n=0; n < nMax; n++){ //forward iteration

		if (cabs2(z)> ER2 || cabs(dz)> 1e60) break; // big values
  			
		dz = derivative_wrt_z(z) * dz; 
		z  = f(z); /* forward iteration : complex quadratic polynomial */ 
	}
  
	if (n==nMax) return iColorOfInterior; // not escaping
	
	
	// escaping and boundary
	cabsz = cabs(z);
	distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
	double g = tanh(distance / PixelWidth);
  
	return 255*g; // iColorOfExterior;
}

• Julia set boundary = DEM/J
• 
• 

logistic map

Logistic map

complex cubic polynomial

wrt z

(%i2) display2d:false;

(%o2) false
(%i3) f(z):= z^3 + c;

(%o3) f(z):=z^3+c
(%i4) diff(f(z), z,1);

(%o4) 3*z^2
(%i5) diff(f(z), c,1);

(%o5) 1

Julia set for f(z) = z^3 +z*(0.1008317508132964*i + 1.004954206930806)

// function 
f(z)= z^3 + c*z 

// first derivativwe wrt z
d = f'(z) = 3*z^2 + c

// iteration:
z_(n+1) = f(z_n)


// initial values:
z = z
d = 1

// Recurrence steps:
d = (3*z^2 + c)*d
z = z^3 + c*z 

quarctic

(%i1) display2d:false;

(%o1) false

(%o4) f(z):=z^4+c
(%i5) diff(f(z), z,1);

(%o5) 4*z^3

FAQ

• why new d must be computed before new z ?

Aplications

• Vector field gradient
• Newton method and recurrence relations
• eDEM = Distance Estimation Method for exterior of the set
  • DEM/M- for Mandelbrot set
  • DEM/J for Julia set
• iDEM = DEM for interior
• stability of periodic point ( multiplier)
• finding critical points
• slope - an algorithm for converting 2D image to 3D

Keep track of derivatives of Z+z  wrt. Z1+z1  (where Z0+z0  is at the critical point, usually 0 ). When the absolute value of the derivative drops below a threshold such as 0.001, classify it as interior and stop iterating.
Keep track of derivatives of Z+z  wrt. pixel coordinates k . As Z  is constant for the whole image, you just need dzdk . An easy way to do this is with dual numbers for automatic numeric differentiation. Set up the pixel coordinates as dual numbers with dual part 1+0i , then transform them to the complex C+c plane of the fractal iterations. At the end you plug the complex derivative into the (directional) distance estimate formula, it is already prescaled by the pixel spacing (this also helps to avoid overflow during iteration). (	Claude Heiland-Allen)

finding equation of normal and tangent to the curve

to find the equation of the normal line to the curve at the given point by Krista King Math:

• Take the derivative of the original curve, and evaluate it at the given point. This is the slope of the tangent line, called m.
• Find the negative reciprocal of m. This is the slope of the normal line, which we’ll call n. So n = −1/m.
• Plug n and the given point into the point-slope formula for the equation of the line, (y−y1)=n(x−x1)
• Simplify the equation by solving for y

Interior distance estimation

The Interior distance estimation is given by :

$distance(c,\sigma M)=d(c,z_{0},p)={\frac {1-\mid {{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}\mid ^{2}}{\mid {{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})+{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0}){\frac {{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})}{1-{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}}}\mid }}$

where

$p$ is the period = length of the periodic orbit

$c$ is the point from parameter plane to be estimated

$f_{c}(z)$ is the complex quadratic polynomial $f_{c}(z)=z^{2}+c$

$f_{c}^{p}(z_{0})$ is the $p$-fold iteration of $f_{c}(z)$

$z_{0}$ is any of the $p$ points that make the periodic orbit ( limit cycle ) $\{z_{0},\dots ,z_{p-1}\}$

4 derivatives of $f_{c}^{p}$ evaluated at $z_{0},c,p$ :

${\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

First partial derivative with respect to z

It must be computed recursively by applying the chain rule :

Starting points : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

First iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

Second iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

General rule for p iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

First partial derivative with respect to c

It must be computed recursively by applying the chain rule :

Starting points : ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

First iteration : ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{0}(z_{0})+1=1\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

Second iteration : ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})+1=2z_{0}^{2}+2c+1\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

General rule for p iteration : ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})+1\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

Second order partial derivative with respect to z

It must be computed :

• together with :$f_{c}^{n}$ and ${\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$
• recursively by applying the chain rule

Starting points : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

First iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0}))^{2}+f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})\}=2\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

Second iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

General rule for p iteration : ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))^{2}+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\}\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

Second order mixed partial derivative

${\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2({\frac {\partial }{\partial {c}}}f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))$

How to compute derivative ?

• wolfram alpha
• Maxima

p-norm of the gradient

The p-norm of the gradient is:

$||g||_{p}=(g_{x}^{p}+g_{y}^{p})^{(1/p)}$

and set p=-2.

See also

• automatic differentiation

derivatives w.r.t. pixel coordinates by Claude

• Pixel = (x + i y); x in [0..w), y in [0..h) (may be fractional if using jitter)
• C = r exp(i pi (Pixel - (w + i h) /2) / w) + C0
• dc/dpixel = r i pi / w exp(i pi (Pixel - (w + i h) /2) / w)
• for z <- z^2 + c:
  • dz/dpixel <- 2 z dz/dpixel + dc/dpixel
  • distance estimate in screen space : de = |z| log |z| / |dz/dpixel|

here:

• pixel_size = |dc/dpixel| (just divide through by dc/dpixel to get dz/dpixel / dc/dpixel = dz/dc, etc)
• dy = (2 * circle_period * Math.PI) / image_size;
• dx = dy * height_ratio;
• height_ratio is usually 1. Its just there to stretch the image on one direction.

derivatives wrt pixel are less likely to overflow, and give screen space de automatically.

Also, now that I checked again the e^(x*dx) term, is actually e^(x*dx + offset). I don't know if this changes the things alot.

 C value = (e^(x*dx + offset)) * (cos(y * dy) + sin(y* dy)i) + C0