• How to draw ...  (example code)

These codes can be written with minor (or no) modifications for java,javascript or c++ to plot:

The Mandelbrot set, the Lorenz attractor, the logistic map/Feigenbaum diagram, the Barnsley fern, a line or 2D surface governed by the wave equation and the Schrödinger equation.

A 'plotPointAt(x,y)' function is required to draw a pixel at (x,y) on the screen at some suitable scale.

N.B. These codes are not optimised for speed, simply for generating a set of points on the screen.

How to draw the Mandelbrot Set

for(i = 0; i < 1000; i++){

  a0 = -2+3*i/1000

  for(j = 0; j < 1000; j++){

    b0 = -1.5+3*j/1000;

    a = a0;

    b = b0;

    for(k=0; k<500; k++){

      as = a0 + a*a-b*b;

      b = b0 + 2*a*b;

      a = as;

      if( a*a + b*b > 3 ){

        colour = factor(k);

        plotPointAt( a0, b0 );

        break;

      }

    }

  }

}

How to draw the Lorenz attractor

// plot 15000 points following the Lorenz attractor

sigma = 6.58;

rho = 50;

beta = 2.67;

x = 1;

y = 1;

z = 1;

dt = 0.002;

for(i = 0; i < 15000; i++){

  xs = x + dt*sigma* (y - x);

  ys = y + dt*x* (rho - z);

  z = z + dt* (x*y - beta*z);

  x = xs;

  y = ys;

  plotPointAt( x, y );

}

How to draw the Logistic Map bifurcation tree / Feigenbaum diagram

// plot 1000 settled points of the logistic map sequence 

// for values of a ranging from 2 to 4

for(i = 0; i < 1000; i++){

  a = 2+i*(4-2)/1000;

  x = 0.1;

  for(j = 0; j < 1000; j++) x=a*x*(1-x);

  for(j = 0; j < 1000; j++){

    x = a*x*(1-x);

    plotPointAt( i, x );

  }

}

How to draw the Barnsley fern

x = 0;

y = 0;

for(i = 0; i < 10000; i++){

  r = random(1);

  if( r<0.02 ) { a=0; b=0; c=0; d=0.16; f=0; }

  else if( r<0.89 ) { a=0.85; b=0.04; c=-0.04; d=0.85; f=1.6; }

  else if( r<0.96 ) { a=0.2; b=-0.26; c=0.23; d=0.22; f=1.6; }

  else { a=-0.15; b=0.28; c=0.26; d=0.24; f=0.44; }

  xs = x;

  x = a*x + b*y;

  y = c*xs + d*y + f;

  plotPointAt( x, y );

}

How to plot the classical wave equation applied to a pulse (in 1D)

// To set up a curved pulse which will then move according to the wave equation:

// Set up p and v as 1D arrays size [100]

for(i = 0; i < 100; i++){

  p[i] = exp(-pow(x-50,2)/100);

  v[i] = 0;

}

// Repeat this calculation at regular time intervals and plot:

for(i = 1; i < 99; i++) v[i]+=p[i-1]-2*p[i]+p[i+1];

for(i = 1; i < 99; i++){

  p[i] += c*v[i];

  plotPointAt(i,p[i]);

} 

How to plot the classical wave equation applied to a bump (in 2D)

// To set up a curved bump on a 2D surface which will then move according to the wave equation:

// Set up p and v as 2D arrays size [100][100]

for(i = 0; i < 100; i++){

  for(j = 0; j < 100; j++){

    p[i][j] = exp((-pow(i-50,2)-pow(j-50,2))/100);

    v[i][j] = 0;

  }

}

// Repeat this calculation at regular time intervals and plot:

for(i = 1; i < 99; i++){

 for(j = 1; j < 99; j++){

  v[i][j]+=-4*p[i][j]+p[i+1][j]+p[i-1][j]+p[i][j+1]+p[i][j-1];

 }

}

for(i = 1; i < 99; i++){

 for(j = 1; j < 99; j++){

  p[i][j] += c*v[i][j];

  plotPointAt(i,j,colour[p[i][j]]);

 }

} 

How to plot the Schrödinger equation for a pulse in 1D

// To set up a curved pulse which will then move according to the Schrödinger equation:

// Set up p_re, p_im, del2Re and del2Im as 1D arrays size [100]

// p (_re and _im) the real and imaginary parts of the wave function, 

// del2Xx the Laplacian operator ∇²(Xx)

//Initial state of wave function   

for(i = 0; i < 100; i++){

  p_re[i] = exp(-pow(x-50,2)/100)*cos(i/3);

  p_im[i] = exp(-pow(x-50,2)/100)*sin(i/3);

}

// Repeat this calculation at regular time intervals and plot:

for(i = 1; i < 99; i++){

 del2Im[i]=p_im[i]*2-p_im[i+1]-p_im[i-1];;

}

for(i = 0; i < 100; i++){

 p_re[i]-=del2Im[i]/2;

}

for(i = 0; i < 100; i++){

 del2Re[i]=p_re[i]*2-p_re[i+1]-p_re[i-1]; 

}

for(i = 0; i < 100; i++){

  p_im[i]+=del2Re[i]/2;

  plotPointAt(i,10+10*p_re[i],red);

  plotPointAt(i,10+10*p_im[i],blue);

} 

How to plot the Schrödinger equation for a bump in 2D

// Initial state of a wave function for a vertically moving object

//  as a plane wave constrained by a circular Gaussian function: 

for(i = 0; i < 100; i++){

 for(j = 0; j < 100; j++){

   p_re[i][j] = exp((-pow(i-50,2)-pow(j-50,2))/100)*cos(j/3);

   p_im[i][j] = exp((-pow(i-50,2)-pow(j-50,2))/100)*sin(j/3);

 }

}

// Repeat this calculation at regular time intervals and plot:

for(i = 1; i < 99; i++){

 for(j = 1; j < 99; j++){

  del2Im[i][j]=p_im[i][j]*4-p_im[i+1][j]-p_im[i-1][j]-p_im[i][j+1]-p_im[i][j-1];

 }

}

for(i = 0; i < 100; i++){

 for(j = 0; j < 100; j++){

  p_re[i][j]-=del2Im[i][j]/2;

 }

}

for(i = 0; i < 100; i++){

 for(j = 0; j < 100; j++){

  del2Re[i][j]=p_re[i][i]*4-p_re[i+1][j]-p_re[i-1][j]-p_re[i][j+1]-p_re[i][j-1];

 }

}

for(i = 0; i < 100; i++){

 for(j = 0; j < 100; j++){

  p_im[i][j]+=del2Re[i][j]/2;

  plotPointAt(i,j,colour[sqrt(pow(p_re[i][j],2)+pow(p_im[i][j]))]);

 }

}