Category Archives: Plotting

How to plot the Julia Set

In the previous post we discussed the Mandelbrot set. There is another beautiful algebraic fractal, called Julia set 𝓙. It is named after the French mathematician Gaston Julia, who has discovered. It is obtained by the same formula as the Mandelbrot Set:

zn+1 = zn2 + c

The difference is that z0 = x + yi is the current point (x, y) in the complex plane and c is a fixed complex number, while for the Mandelbrot set z0 = 0 and c = x + yi.

Both sets are closely related. There are a lot of different Julia sets, corresponding to different values of c. Their shapes and properties depend on where the point is located, related to the Mandelbrot set 𝓜. Inside 𝓜, the corresponding Julia set 𝓙 is connected. Outside 𝓜, it falls apart into a set of isolated points, called Fatou set 𝓕. Further away, it turns into Fatou dust – infinite number of points with zero area, similarly to the Cantor dust.

For simple plotting, you can use the following code:

"Julia set
'Define the function
JuliaSet(z; c) = $Repeat{z = z^2 + c @ i = 1 : n}
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 600','PlotHeight = 400
#show
'Plot for'c = -0.4 + 0.59i'and'n = 200'iterations
$Map{abs(JuliaSet(x + 1i*y; c))) @ x = -1.5 : 1.5 & y = -1 : 1}

When you run it inside Calcpad, you will get the following result:

Plotting Julia set with Calcpad

To obtain better and colorful images, we will use different approach this time. We will modify the JuliaSet function in a way to return directly the number of iterations for which the result overpasses 2. The source code is provided bellow:

"Julia set
J(z; c; n) = abs($Repeat{z = if(abs(z) < 2; z^2 + c; 4) @ i = 1 : n})
JuliaSet(z; c) = $Find{J(z; c; x) - 2 @ x = 1 : n}
#hide
PlotStep = 1','PlotWidth = 600','PlotHeight = 400','Precision = 10^-6
#show
'Plot for'c = -0.4 + 0.6i'and'n = 100'iterations
$Map{JuliaSet(x + 1i*y; c)) @ x = -1.5 : 1.5 & y = -1 : 1}

It runs a little bit slower, but the result is worth the wait:

Colorful plot of Julia set

You can download Calcpad for free and experiment by yourself with different values of c. For example, the lightning image bellow is a special case of Julia set for c = i and is called “dendrite” fractal.

Dendrite fractal image, obtained for c = i

Plotting the Mandelbrot Set

Calcpad has an interesting and undocumented feature, that I am going to reveal in this post. You can use it to quickly plot the Mandelbrot set. This is a set of complex numbers c, for which the iterative equation zn+1 = zn2 + c does not go to infinity. The most beautiful part is that it shows fractal behavior at its boundaries. If you plot it with the appropriate colors, you can create stunning images and animations. I really like those on the Maths Town’s YouTube channel:

Mandelbrot fractal zoom videos on YouTube

Usually, such videos are created by complicated software, specially developed for that purpose. However, you can also plot the Mandelbrot set with Calcpad and a few lines of code. There are two ways to do that. The first one is to use the Calcpad capabilities for complex arithmetic, as follows:

"Mandelbrot set
'Define the function
MandelbrotSet(z; c) = $Repeat{z = z^2 + c @ i = 1 : n}
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 500','PlotHeight = 500
#show
'Plot for'n = 50'iterations
$Map{abs(MandelbrotSet(0; x + 1i*y)) @ x = -1.5 : 0.5 & y = -1 : 1}

Open Calcpad, switch to “complex”mode, paste the above code inside, and run the calculations. You will get the following result:

Mandelbrot set plot with Calcpad complex arithmetic

It is not too bad, but still not very impressive. At least, the exterior is missing any colors. Then, how we can do the coloring? Theoretically, it is satisfied for all numbers in the set that |zn| < 2 for any n. Outside, the process diverges and the modulus goes to infinity, as n increases. You can use that in the following way:
1. For each point, calculate the number of iterations n, where |zn| gets larger than 2.
2. Map these values to colors.
3. Plot the colors at the respective coordinates.

If you play with the number of iterations n in the above example, you will see that the plot changes. The shaded area is where |zn| converges for the respective n.

Mandelbrot set, plotted with different values of n

If you combine the above plots and apply some colors, you can get really beautiful artistic image like the one bellow:

Mandelbrot set with Calcpad + GIMP

The only problem with this method is that it is slow and elaborate. That is why, I included a special function Mandelbrot(x; y) in Calcpad. It calculates the number of iterations, needed for the modulus to overpass 2. You can plot this function directly with the preferred color scale:

Mandelbrot function plot
Mandelbrot function partial plot

To do that, I used even shorter piece of code:

"Mandelbrot set
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 500','PlotHeight = 500
#show
'Plot for'n = 50'iterations
$Map{Mandelbrot(x; y) @ x = -1.5 : 0.5 & y = -1 : 1}
'Plot for'n = 70'iterations
$Map{Mandelbrot(x; y) @ x = -0.59 : -0.58 & y = 0.55 : 0.56}

Here you can find the C# code of the Mandelbrot function that works inside Calcpad:

        private static readonly double Log2Inv = 1 / Math.Log(2);       
        //Calculates the number of iterations for which z = z^2 + c satisfies |z| > 2
        //Then transforms it to a smooth log scale and returns the result
        protected static double MandelbrotSet(double x, double y)
        {
            //Checks if the point is inside the set and directly returns NaN
            if (x > -1.25 && x < 0.375)
            {
                if (x < -0.75)
                {
                    if (y > -0.25 && y < 0.25)
                    {
                        double x1 = x + 1,
                            y2 = y * y,
                            x2 = x1 * x1;
                        if (x2 + y2 <= 0.0625)
                            return double.NaN;
                    }
                }
                else if (y > -0.65 && y < 0.65)
                {
                    double x1 = x - 0.25,
                        y2 = y * y,
                        q = x1 * x1 + y2;
                    if (q * (q + x1) <= 0.25 * y2)
                        return double.NaN;
                }
            }
            //For all other points performs detailed calculations
            double re = x, im = y;
            for (int i = 1; i <= 1000; ++i)
            {
                double reSq = re * re, 
                       imSq = im * im,
                       sumSq = reSq + imSq;
                //To avoid the sqrt function, the check |z| > 2 is replaced by z^2 > 4
                if (sumSq > 4)
                {
                    var logZn = Math.Log(sumSq) / 2;
                    var nu = Math.Log(logZn * Log2Inv) * Log2Inv;
                    return (1.01 - Math.Pow(i - nu, 0.001)) * 1000;
                }
                //Calculates z = z^2 + c
                im = 2 * re * im + y;
                re = reSq - imSq + x;
            }
            return double.NaN;
        }
    }

The complete source code is available on GitHub.

Even more 3D effects

In the last version I improved the 3D plot shading and added some nice 3D effects. Bellow is the plot of the function

f(x; y) = cos(x) + cos(y) + cos(x·y/3)/2 + cos(√(x2 + y2)·2)/2

with the “Rainbow” color scale and the “Shadows” option on:

Specular

What is new is the specular light effect (the little shiny reflections). I also replaced the Gouraud shading procedure with the Phong one. It improved the image quality  without any noticeable impact on the speed (surprisingly). Since the light is distant and the view is isometric, I use the Blinn–Phong shading model. Instead of calculating the reflection vector for each vertex, it calculates the halfway vector between the viewer and light-source. It is performed once for the whole model, and that is why it is faster.

In general, there is additional overhead for generating the Html output, converting the image to base64 and rendering the Html into the “Output” window. However, as long as it takes much less than a second for all of it, the speed is not of a big concern.

Math Art with Calcpad

Most people consider math boring, but sometimes it can be really beautiful. Recently I found some unusual application of Calcpad – to draw nice pictures using math formulas. In general, you can plot functions of two variables using the $Map command:

$Map{f(x;y) @ x = a : b & y = c : d}

I notices that you can get nice effects using the following types of functions:

f(x;y) = abs(cos(p1(x;y)) + cos(p2(x;y)) ...)

It creates families of intersecting curves that from canyons with hills in between. If you select “none” for the color scale, it creates nice and smooth pictures with shadows only. Bellow are some examples:

circles-and-hyperbolas

f(x;y) = abs(cos(x^2/10) + cos(y^2/10))

unsymmetric-hyperbolas

f(x;y) = abs(cos(x + y/3) + cos(x*y/7))

zig-zag-tiles

f(x;y) = abs(cos(x – sin(y)) + cos(y – sin(x)))

concentric-waves

f(x;y) = abs(cos(r(x;y)*e^(ξ*r(x;y)))*e^(-1.5*ξ*r(x;y)))

intersecting-parabolas

f(x;y) = abs(cos(r(x;y)) + cos(y))

intersecting-parabolas-x

f(x;y) = abs(cos(x) + cos((y/4)^2))

puzzle

f(x;y) = abs(cos(2*r(x;y)) + cos(x) + cos(y))

buttons

f(x;y) = abs(cos(2*r(cos(x/2);cos(y/2))))

pins

f(x;y) = abs(r(cos(x/2);cos(y/2)))

In the above equations, r(x;y) = sqr(x^2 + y^2).

Download the latest Calcpad 3.2 and try to create your own pictures. Enjoy!

Plotting with Calcpad 3.2

The new and enhanced version of Calcpad 3.2 was released. Besides other improvements, we added some nice visual effects for 2D plotting.

In general, you can use Calcpad to plot functions of two parameters f(x; y), by the $Map command:

$Map{f(x;y) @ x = a : b & y = c : d}

The result is displayed as a 2D contour map. For example, let us take the function:

f(x; y) = cos(x) + cos(y) + cos(π/2*sqr(x^2 + y^2))

The $Map command will produce the following image:

The left one is without shadow effects. Different values of the function are represented by the respective colors. Blues are the lowest and reds are the highest. However, you need to have quite a rich imagination to figure out how the surface actually looks. That is why we added a 3D shadows effects. On the right side, you can see the same image, but with shadows. The light here comes from North-West direction. Now it looks much more natural and it is easier for the human eye to perceive the shape.

The shading is calculated using the dot product of the normal vector to the surface by the vector of the light.

There are also other plotting options available:

Toolbar

You can select among different color scales:

If you select “smooth” (gradient) scale, you will get the following image:

Gradient map with shadows

Instead of having contours, now the colors flow smoothly from one to another. That is because the color of each pixel is calculated individually based on its value.

Download the latest version of Calcpad 3.2 and try by yourself:

https://calcpad.eu/download/calcpad-setup-en-x64.zip