Julia Sets

Julia Sets are another set of images that can be created by iterating a dynamic equation. If you want all of the ugly mathematics, see the Wikipedia entry on Julia Set. If you take the time to understand all of it, you are a better person than me. For a better explanation, I suggest seeing Julia and Mandelbrot Sets, and Understanding Julia and Mandelbrot Sets.

That’s all fine and good, but how do we create a Julia Set mathematically? Certainly not by attempting to follow the Wikipedia entry! For a Julia Set that corresponds to the Mandelbrot set, we use the same equation, but generate the plot for a different plane. Remember the equation:

zn+1 = znp + c

For the Mandelbrot Set, z0 = 0 + 0i and p = 2. We set c to be each point in the section of the complex plane that we were looking at. For the Julia Set, we set c to a single value in the complex plane, then set z0 to each point in the section of the complex plane we are looking at, and iterate. For every point, c, inside the Mandelbrot Set, there is a connected Julia Set; for each point outside of the Mandelbrot Set, the Julia Set is disconnected.

The most interesting images are generated from points outside the Mandelbrot Set or just inside it. This image is generated for a point near the centre of the secondary bulb:

MJulia1

and this is image is generated for a point inside one of the small bulbs off the main bulb:

MJulia2

This image is generated for a point just outside the Mandelbrot set and slightly about the x-axis:

MJulia3

Here are a couple of images created from the tendrils of the Mandelbrot Set:

MJulia4

MJulia5

and finally, an image from a short distance outside the Mandelbrot Set:

MJulia6

Each value, c, produces a different image. To generate an image for a point in the Mandelbrot Set using the ChaosExplorer program, place the mouse cursor over the point, click the right mouse button, and select the Julia Set menu item. The image is generated in a panel in a new notebook tab. You can also zoom in on the image by selecting an area to zoom, and then selecting the Draw From Selection menu item. The source code and a wiki are available on GitHub.

This is the fragment shader for generating these Julia Sets:

    std::string fragmentSource =
        "#version 330 core\n"
        "uniform vec2 c;"
        "uniform vec2 p;"
        "uniform vec2 ul;"
        "uniform vec2 lr;"
        "uniform vec2 viewDimensions;"
        "uniform vec4 color[50];"
        "out vec4 OutColor;"
        "vec2 iPower(vec2 vec, vec2 p)"
        "{"
        "    float r = sqrt(vec.x * vec.x + vec.y * vec.y);"
        "    if(r == 0.0f) {"
        "        return vec2(0.0f, 0.0f);"
        "    }"
        "    float theta = vec.y != 0 ? atan(vec.y, vec.x) : (vec.x < 0 ? 3.14159265f : 0.0f);"
        "    float imt = -p.y * theta;"
        "    float rpowr = pow(r, p.x);"
        "    float angle = p.x * theta + p.y * log(r);"
        "    vec2 powr;"
        "    powr.x = rpowr * exp(imt) * cos(angle);"
        "    powr.y = rpowr * exp(imt) * sin(angle);"
        "    return powr;"
        "}"
        ""
        "void main()"
        "{"
        "    float x = ul.x + (lr.x - ul.x) * gl_FragCoord.x / (viewDimensions.x - 1);"
        "    float y = lr.y + (ul.y - lr.y) * gl_FragCoord.y / (viewDimensions.y - 1);"
        "    vec2 z = vec2(x, y);"
        "    int i = 0;"
        "    while(z.x * z.x + z.y * z.y < 4.0f && i < 200) {"
        "        z = iPower(z, p) + c;"
        "        ++i;"
        "    }"
        "	 OutColor = i == 200 ? vec4(0.0f, 0.0f, 0.0f, 1.0f) : "
        "        color[i%50 - 1];"
        "}";

If you compare this to the fragment shader for the Multibrot Sets included in the post Mandelbrot Set, you will see few changes.

Summary

This post discussed Julia Sets and showed images of some Julia Sets generated from several points inside and outside the Mandelbrot Set. The fragment shader that is used to generate the images was also shown.

Advertisements

Mandelbrot Set

Complex Numbers and the Complex Plane

Before we discuss and generate the Mandelbrot set, we need to look at complex numbers, the complex plane, and Euler’s formula. If you are familiar with these concepts, then you can skip to the next section in this post. Or, if you can’t stand the math, just go look at the pretty pictures. ūüôā

A complex number is of the form a + bi, where a and b are real numbers and i is defined such that i2 = -1. Note that i is defined only in terms of what the value of i2 is, and not by itself. You might quickly think that i = sqrt(-1), but this is only one of the roots of the equation. In fact, as my university algebra professor showed, by defining i this way, it is possible to prove that -1 = +1. That would certainly negate all of mathematics (and physics) as we know it.

Complex numbers are required to solve equations for which there is no real solution. Take for example:

(x + 1)2 = -9

There are no real number values for x that satisfy that equation. However, as shown here, there are two complex number values that do satisfy the equation:

x = -1 + 3i, and

x = -1 -3i.

A complex number can be viewed as a point or a position vector plotted in a two-dimensional Cartesian coordinate system. This is called the complex plane. Complex numbers are normally plotted such that the real portion of a complex number x + yi is plotted with the x value along the horizontal axis, and the y value plotted along the vertical axis.

Now that we have plotted a complex number in the two-dimensional Cartesian coordinate system, it is possible to specify the complex number using polar coordinates, and finally, using Euler’s formula.

See the Wikipedia page on Complex numbers for¬†more information. Using Euler’s formula, it is even possible to calculate the complex power of a complex number. The Wolfram Mathworld website contains a web page that shows the formula for calculating these values. And finally, the Abecedarical Systems website contains a web page with¬†the listing of a program that calculates the complex power of a complex number.

Now, why do we need all of this mathematics? C++ contains a templated complex number class, but GLSL does not support complex numbers. Therefore, to calculate the complex power of a complex number using GLSL, it will be necessary to use Euler’s formula. You will see this later when we look at the GLSL code for the fragment shader.

Mandelbrot Set

Now that we have covered some rather heavy mathematics, we can look at the Mandelbrot set. Assuming z0 = 0 and c is a point in the complex plane, then the Mandelbrot set is the set of all points c in the complex plane where the value zn+1 does not tend toward infinity as n approaches infinity in the following equation:

zn+1 = zn2 + c.

If we attempt to iterate this formula an infinite number of times for each point in the complex plane, the calculations will of course take an infinite amount of time. There are two pieces of information that we an use to limit the number of calculations:

  1. For all values of zn where |zn| ‚Č• 2, the value approaches infinity as n approaches infinity.
  2. As n increases, the set of points approaches the Mandelbrot set quite quickly.¬†There is little difference between the points defined by z200 and z‚ąě.

The program, ChaosExplorer, uses these two facts to limit the number of iterations that are required for each point. ChaosExplorer is a C++ program that uses wxWidgets for the windowing system, and GLSL to calculate and display the points that are inside and outside the Mandelbrot set. The majority of the code will not be shown or discussed in this post because similar or identical code has been covered in earlier posts. The program code is available on GitHub.

When ChaosExplorer starts, it displays the Mandelbrot set, shown below, in a tabbed window in the main window of the program.

Mandelbrot

 

This image shows the complex plane for -2.5 < x < 1.5 and -2.0 < y < 2.0. The portion in black is the Mandelbrot set. The colours change as the number of iterations (n) required for zn to exceed the value 2 increases. There are 50 colours used, with the full set repeated 4 times for a total of 200 iterations.

You can select an area to zoom in on by pressing the left mouse button and dragging towards the right and down. As you select an area, the area is highlighted by a white line around the selection. When you have chosen the area you want to select, release the left mouse button. A selection is shown in the image above. Now to zoom into that area, click the right mouse button anywhere in the image to display the context menu, and select the Draw From Selection menu item. The resulting display is shown in a second tabbed window. Here is the result for the region selected in the image above:

Mandelbrot2

Note the smaller version of the main Mandelbrot set. You can continue to zoom in on the tendrils, such as the next image, to see ever smaller versions of this set. This shows that the Mandelbrot set is a fractal.

You can also animate the zoom by placing the mouse cursor at the point you want to zoom, clicking the right mouse button, and selecting Animate Magnification from the context menu.

Here are some additional zoomed images:

Mandlebrot4

Mandelbrot3

Mandelbrot5

and a video of a zooming image:

To redisplay the full Mandelbrot set, select the Power=2 menu item in the MultibrotPower submenu.

Above, I stated that points selected by each iteration quickly approaches the Mandelbrot set as the number of iterations increases. You can see this by selecting the Animate Iterations menu item in the context menu. The image will change for each iteration, showing the pixels that are excluded by each iteration. You can run this animation at any zoom level.

I leave it to you as an exercise to investigate the various regions of the image at varying magnifications.

So, how are the images calculated and drawn? The calculations are done in the fragment shader:

    std::string fragmentSource =
        "#version 430 core\n"
        "uniform vec2 z0;"
        "uniform vec2 p;"
        "uniform vec2 ul;"
        "uniform vec2 lr;"
        "uniform vec2 viewDimensions;"
        "uniform int maxIterations;"
        "uniform vec4 color[50];"
        "out vec4 OutColor;"
        "vec2 iPower(vec2 vec, vec2 p)"
        "{"
        "    float r = sqrt(vec.x * vec.x + vec.y * vec.y);"
        "    if(r == 0.0f) {"
        "        return vec2(0.0f, 0.0f);"
        "    }"
        "    float theta = vec.y != 0 ? atan(vec.y, vec.x) : (vec.x < 0 ? 3.14159265f : 0.0f);"
        "    float imt = -p.y * theta;"
        "    float rpowr = pow(r, p.x);"
        "    float angle = p.x * theta + p.y * log(r);"
        "    vec2 powr;"
        "    powr.x = rpowr * exp(imt) * cos(angle);"
        "    powr.y = rpowr * exp(imt) * sin(angle);"
        "    return powr;"
        "}"
        ""
        "void main()"
        "{"
        "    vec2 z = z0;"
        "    float x = ul.x + (lr.x - ul.x) * gl_FragCoord.x / (viewDimensions.x - 1);"
        "    float y = lr.y + (ul.y - lr.y) * gl_FragCoord.y / (viewDimensions.y - 1);"
        "    vec2 c = vec2(x, y);"
        "    int i = 0;"
        "    while(z.x * z.x + z.y * z.y < 4.0f && i < maxIterations) {"
        "        z = iPower(z, p) + c;"
        "        ++i;"
        "    }"
        "	 OutColor = i == maxIterations ? vec4(0.0f, 0.0f, 0.0f, 1.0f) : "
        "        color[i%50 - 1];"
        "}";

GLSL does not contain a complex data type, so complex values are passed in as vec2 uniform values. The iPower function calculates the complex power of a complex number using Euler’s formula as mentioned in the first section of this post. For integer real powers of complex numbers, the mathematics could be simplified considerably, but the full calculation in iPower is included because it will be used in other functionality in this program.

Finally, the main function performs repeated iterations of the formula zn+1 = znp + c until the magnitude of zn+1 exceeds the escape value of 2. c is the complex number that specifies the coordinates of the point in the complex plane. The color for the fragment is then determined from the number of iterations (the value of i in the code). This video shows the full Mandelbrot set as the number of iterations increases. The number of iterations is shown in the status bar at the bottom of the window. Once the number of iterations exceeds about 40, the only place you will see any changes at this zoom level is at the junction between the main bulb and the bulb to its left.

Advantages and Disadvantages of Performing the Calculations in the Fragment Shader

Advantages

The main advantage is the speed at which the image is calculated and drawn. Using modern GPUs with hundreds or thousands of cores, the image can be updated at the rate of many frames per second as seen when animating the magnification. If the image was calculated in the CPU, the number of frames per second that could be generated is much less (around 1 frame per second) even if the code were multithreaded.

Disadvantages

  1. GLSL does not contain a complex data type, so Euler’s formula must be used to calculate a¬†complex power of a complex number.
  2. Since version 4.0.0 of OpenGL, GLSL supports the double data type, but trigonometric functions are supplied only for floating point values. This means that the complex numbers can be specified only as two floating point numbers and not as two doubles. As the zoom factor of the image increases. the limitations on the values of floating point numbers become a factor, and pixelation of the image occurs much sooner than would be the case if doubles could have been used. Note the pixelation in the following image:

pixelated

Summary

This post discussed the Mandelbrot set and the ChaosExplorer program. Several images of the Mandelbrot set are shown, as are examples of the Animate Iterations and Animate Magnification functions of the program.

The Mandelbrot set is a specialization of the formula zn+1 = znp + c where p = 2. The next post will look at the more general Multibrot set, wherein p takes on other values, including complex values.

More information about ChaosExplorer and the Mandelbrot Set are included in the wiki for the ChaosExplorer program on GitHub.

Chaotic Systems and Fractals

Chaos Theory

Chaos theory is the study of dynamical systems that are highly sensitive to initial conditions. So, a chaotic system is one in which small changes in initial conditions result in large variations in outcomes. Weather is one example of a chaotic system. At the 139th meeting of the American Association for the Advancement of Science in 1972, Philip Merilees presented a talk titled Does the flap of a butterfly‚Äôs wings in Brazil set off a tornado in Texas? The hypothesis is that the beating of a butterfly’s wings could create tiny changes in atmospheric conditions that may alter the path of a tornado, or accelerate, delay or prevent a tornado elsewhere.

There are many more examples of chaotic systems:

  • Lorenz Attractor;
  • Weather Forecasting;
  • Double Pendulum;
  • Three Body Gravitational Problem;
  • Randomly Oscillating Magnetic Pendulum;
  • Game of Life;
  • Billiards on an Oval Table (Bunimavich Stadium);
  • Dripping Honey;
  • A Dripping Faucet;
  • Smale Horseshoe;
  • Traffic Flow;
  • Organizational Behaviour;
  • Shifts in Public Opinion; and,
  • Epidemics.

More information on chaotic systems and the butterfly effect are available on the following web pages and lectures:

Fractals

Many chaotic systems display fractal behaviour. A fractal is an infinitely complex pattern that is self-similar at different scales of magnification. A fractal is created by repeating a simple process over and over in an ongoing feedback loop (from What is Chaos Theory?). There are many physical and mathematical examples of fractals; physical examples include:

  • Snowflakes;
  • Ferns;
  • Ice crystals on a window;
  • Romanesco Broccoli;
  • Pine Cone Seeds;
  • Tree Branches;
  • Mountains;
  • Rivers; and,
  • Leaves

to mention but a few.

There are also many mathematical examples of fractals, a few of which we will explore in a program called ChaosExplorer and the next several posts. Note that many or all of the examples have been programmed elsewhere; the difference between ChaosExplorer and the other programs is that ChaosExplorer is an OpenGL program that performs the mathematics in the fragment shader in the GPU; most other programs perform the mathematics in the CPU.

There are both advantages and disadvantages to doing the calculations in the GPU. These will be discussed in the post on the Mandelbrot Set as it is easier to show examples than to simply list them.