{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "etfCiiuVyFrY" }, "source": [ "\n", "\n", "\n", "# Fractals\n", "\n", "---\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " \n", " Dr Jon Shiach

\n", " Senior Lecturer
\n", " Department of Computing and Mathematics
\n", " Manchester Metropolitan University
\n", " Email: j.shiach@mmu.ac.uk\n", " \n", " 		\n", " \n", " \n", " \n", " Malcolm Connolly

\n", " Statistics Tutor
\n", " Department of Computing and Mathematics
\n", " Manchester Metropolitan University
\n", " Email: m.connolly@mmu.ac.uk\n", " \n", " 		\n", " \n", " \n", " \n", " Dr Stephen Lynch

\n", " Reader
\n", " Department of Computing and Mathematics
\n", " Manchester Metropolitan University
\n", " Email: s.lynch@mmu.ac.uk\n", " \n", " 		\n", " \n", " \n", " \n", " Dr Killian O'Brien

\n", " Senior Lecturer
\n", " Department of Computing and Mathematics
\n", " Manchester Metropolitan University
\n", " Email: k.m.obrien@mmu.ac.uk\n", " \n", " 		\n", " \n", "
\n", "\n", "---\n", "\n", "\n", "## Introduction\n", "\n", "[**Fractals**](https://en.wikipedia.org/wiki/Fractal) are complex and beautiful shapes that are generated using a simple set of mathematical rules. Fractals often exhibit a complicated structure which is repeated at different scales throughout the fractal, this means that if we zoom in on a portion of the fractal it will resemble the fractal itself. Mathematicians refer to this property as **self similarity** and can be seen in nature, some examples of which are shown below.\n", "\n", "| | | | |\n", "|:--:|:--:|:--:|:--:|\n", "| Cactus | Electricity | Snowflake | Lungs |\n", "\n", "This Jupyter notebook will introduce you to fractals and how we can draw fractals using [Python](https://www.python.org/). To be able to follow the content of this notebook you just need an understanding of basic trigonometry, geometry and complex numbers (which are covered here). You will be asked to enter Python programs into code cells and run the programs. Typing the programs into the code cells will help you to understand what the program is doing, but if this gets tedious you can also copy and paste the programs into the code cells.\n", "\n", "### Jupyter notebooks\n", "\n", "Jupyter notebooks are documents that combine text and Python code which allow readable documents such as this one to be created that contain executable code used to perform calculations. To run code in a notebook simply click on the code and then on the run button, or alternatively, press the **ctrl + enter** keys at the same time. Since we will be using commands to perform calculations and produce plots of fractals we need to import some commands to help us to do this. Run the code below to import the commands." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "4TjsHTX5yFrc" }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib import cm, colors\n", "plt.rcParams['figure.figsize'] = [10, 8]" ] }, { "cell_type": "markdown", "metadata": { "id": "RS56Fc2-yFrd" }, "source": [ "These commands import the Numpy (short for *Numerical Python*) and matplotlib libraries which enable us to perform numerical calculations and produce plots of our fractals.\n", "\n", "Now that we have imported these libraries we can perform some calculations. The Python code below calculates the length of the side of the right-angled triangle opposite the angle $\\theta = 0.6435$ and with a hypotenuse of length 5 and prints the result. Note how the equation $\\textsf{opposite}=5\\sin(0.6435)$ is entered in a similar way to how we write it on a piece of paper. Can you add a couple of lines of code to calculate the length of the adjacent side as well?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "wC0__JTsyFrd" }, "outputs": [], "source": [ "opposite = 5 * np.sin(0.6435)\n", "print(opposite)" ] }, { "cell_type": "markdown", "metadata": { "id": "dJRWLKBDyFre" }, "source": [ "---\n", "## The Koch Curve\n", "\n", "The [**Koch Curve**](https://en.wikipedia.org/wiki/Koch_snowflake) first discovered by Swedish mathematician Helge von Koch is one of the first examples of a fractal. To draw a Koch curve we begin with a straight line of length $L$ (the initial starting point is known as stage 0). We then remove the middle third of the line and replace it with two lines of length $\\frac{L}{3}$ which form two sides of an equilateral triangle. So after the first time we do this (stage 1) we have four lines each of length $\\frac{L}{3}$. We then repeat the process for each of the new lines.\n", "\n", "\n", "\n", "The Python program below defines a function called `kochCurve` which uses inputs of the starting co-ordinates `x1` and `y1`, the length of the line `L`, the angle of the line `angle` and the number of stages `n`. Run the code cell to see the result." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "OAt0n8dmyFre", "scrolled": false }, "outputs": [], "source": [ "# This is a code cell which we can write Python programs.\n", "# Enter the program here and then run the code cell by clicking on 'Run' or press ctrl and enter.\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "def kochCurve(x, y, L, angle, n):\n", " if n == 0:\n", " # Add the end point co-ordinates to x and y\n", " x.append(x[-1] + L * np.cos(np.radians(angle)))\n", " y.append(y[-1] + L * np.sin(np.radians(angle)))\n", " if n > 0:\n", " L = L / 3;\n", " kochCurve(x, y, L, angle, n - 1) # add line\n", " kochCurve(x, y, L, angle + 60, n - 1) # turn 60 degrees to the left and add line\n", " kochCurve(x, y, L, angle - 60, n - 1) # turn 60 degrees to the right and add line\n", " kochCurve(x, y, L, angle, n - 1) # add line\n", "\n", "\n", "# Generate Koch curve\n", "initialLength = 1;\n", "initialAngle = 0;\n", "numberStages = 2;\n", "x, y = [0], [0]\n", "kochCurve(x, y, initialLength, initialAngle, numberStages)\n", "\n", "# Plot curve\n", "fig, ax = plt.subplots()\n", "plt.plot(x, y, color=\"black\", linewidth=0.5)\n", "ax.axis(\"off\")\n", "ax.set_aspect(\"equal\")" ] }, { "cell_type": "markdown", "metadata": { "id": "9Mg3pZqjyFrf" }, "source": [ "If everything has gone well you should see the Koch curve drawn to stage 2. Can you edit the program to draw the Koch curve up to stage 4?" ] }, { "cell_type": "markdown", "metadata": { "id": "TkvuCgbPyFrf" }, "source": [ "### The Koch snowflake\n", "\n", "The **Koch snowflake** is a fractal that is constructed using a Koch curve to each side of an equilateral triangle.\n", "\n", "\n", "\n", "The Python code below uses our `kochCurve` function that we defined above to draw a Koch snowflake." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "EFUm7riQyFrf", "scrolled": false }, "outputs": [], "source": [ "# Generate the Koch snowflake\n", "initialLength = 1\n", "initialAngle = 0\n", "numberStages = 6\n", "x, y = [0], [0]\n", "kochCurve(x, y, initialLength, initialAngle, numberStages)\n", "kochCurve(x, y, initialLength, initialAngle + 240, numberStages)\n", "kochCurve(x, y, initialLength, initialAngle + 120, numberStages)\n", "\n", "# Plot curve\n", "fig, ax = plt.subplots()\n", "plt.plot(x, y, color=\"black\", linewidth=0.5)\n", "plt.axis(\"equal\")\n", "plt.axis(\"off\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "fAgyYklYyFrg" }, "source": [ "---\n", "### Activity 1: The Koch square\n", "\n", "A variation on the Koch snowflake is the **Koch Square** which is where each line of a square has the middle third removed and replaced by three lines which form three sides of a square with side length one third of that of the line we started with. Amend the program below by replacing the `????` with commands so that is can draw the Koch square." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "E_aECKk-yFrg", "scrolled": false }, "outputs": [], "source": [ "def kochSquare(x, y, L, angle, n):\n", " if n == 0:\n", " # Add the end point co-ordinates to x and y\n", " x.append(x[-1] + L * np.cos(np.radians(angle)))\n", " y.append(y[-1] + L * np.sin(np.radians(angle)))\n", " if n > 0:\n", " L = L / 3; # divide length of line by 3\n", " kochSquare(x, y, L, angle, n - 1) # add first line\n", " # kochSquare(x, y, L, ????, n - 1) # add second line\n", " # kochSquare(x, y, L, ????, n - 1) # add third line\n", " # kochSquare(x, y, L, ????, n - 1) # add fourth line\n", "\n", "\n", "# Generate the Koch snowflake\n", "initialLength = 1\n", "initialAngle = 0\n", "numberStages = 6\n", "x, y = [0], [0]\n", "kochSquare(x, y, initialLength, initialAngle, numberStages) # add top kochSquare curve\n", "# kochSquare(x, y, initialLength, ????, numberStages) # add right-hand side kochSquare curve\n", "# kochSquare(x, y, initialLength, ????, numberStages) # add bottom kochSquare curve\n", "# kochSquare(x, y, initialLength, ????, numberStages) # add left-hand side kochSquare curve\n", "\n", "# Plot curve\n", "fig, ax = plt.subplots()\n", "plt.plot(x, y, color=\"black\", linewidth=0.5)\n", "plt.axis(\"equal\")\n", "plt.axis(\"off\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "DPonpu-JyFrh" }, "source": [ "---\n", "## The Sierpinski Triangle\n", "\n", "The [**Sierpinski Triangle**](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_triangle) was discovered by Polish mathematician Wacław Sierpiński and is drawn by starting with a filled equilateral triangle at stage 0. We then remove the inverted equilateral triangle with the vertices at the midpoints along the edges of the triangle. We repeat this for each new triangle.\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "If the initial triangle has side lengths 1 then the co-ordinates of the three vertices anti-clockwise starting at the bottom-left vertex are $(0, 0)$, $(1, 0)$ and $( \\frac{1}{2}, \\frac{\\sqrt{3}}{2} )$. The vertices of the inverted triangle are labeled $(x_1, y_1)$, $(x_2, y_2)$ and $(x_3, y_3)$ going anti-clockwise starting at the bottom vertex as shown in the diagram on the right. To calculate these co-ordinates we use the co-ordinates of the bottom-left vertex $(x, y)$ and the side length $L$ to give\n", "\n", "\\begin{align*}\n", " (x_1, y_1) &= \\left(x + \\frac{1}{2}L, y\\right), &\n", " (x_2, y_2) &= \\left(x + \\frac{3}{4}L, y + \\frac{\\sqrt{3}}{4}L\\right), &\n", " (x_3, y_3) &= \\left(x + \\frac{1}{4}L, y + \\frac{\\sqrt{3}}{4}L\\right).\n", "\\end{align*}\n", "\n", "The three smaller triangles are defined by the co-ordinates of the bottom-left vertex which are $(x, y)$, $(x_1, y_1)$ and $(x_3, y_3)$ and the side length one-third of the current triangle. The Python program below plots the Sierpinski triangle to stage 2. Edit the program to generate the Sierpinski triangle using more stages." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "KCjSzMK_yFrh" }, "outputs": [], "source": [ "def sierpinskiTriangle(x, y, L, n):\n", " # Calculate co-ordinates of the inverted triangle\n", " x1 = x + 1/2 * L\n", " x2 = x + 3/4 * L\n", " x3 = x + 1/4 * L\n", " y1 = y\n", " y2 = y + np.sqrt(3) / 4 * L\n", " y3 = y2\n", "\n", " # Plot inverted triangle\n", " plt.fill([x1, x2, x3], [y1, y2, y3], facecolor=\"white\")\n", "\n", " # Repeat for each sub-triangle\n", " if n > 1:\n", " L = L / 2\n", " sierpinskiTriangle(x, y, L, n - 1) # bottom left triangle\n", " sierpinskiTriangle(x1, y1, L, n - 1) # bottom right triangle\n", " sierpinskiTriangle(x3, y3, L, n - 1) # top triangle\n", "\n", "\n", "# Plot stage 0 triangle\n", "initialLength = 1\n", "fig, ax = plt.subplots()\n", "plt.fill([0, initialLength, initialLength / 2], [0, 0, np.sqrt(3) / 2], facecolor=\"blue\")\n", "plt.axis(\"equal\")\n", "plt.axis(\"off\")\n", "\n", "# Generate Sierpinsky triangle\n", "numberStages = 2\n", "sierpinskiTriangle(0, 0, initialLength, 2)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "C-Fq2raqyFrh" }, "source": [ "---\n", "### Activity 2: The Sierpinski carpet\n", "\n", "The [**Sierpinski carpet**](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_carpet) is similar to the Sierpinski triangle but drawn by starting with a filled square at stage 0. We divide the square into 9 squares and remove the centre square. We repeat this step for each of the 8 filled squares of side lengths one third of the current square that remain. Amend the program below so that it can draw the Sierpinski carpet.\n", "\n", "\n", "\n", "Amend the program below by replacing `????` with appropriate commands so that it can draw the Sierpinsky carpet." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "SklGMMKDyFri", "scrolled": false }, "outputs": [], "source": [ "def sierpinskiCarpet(x, y, L, n):\n", " # Calculate co-ordinates of the inverted square\n", " x1 = x + 1/3 * L # x co-ordinate of lower left corner of middle square\n", " y1 = y + 1/3 * L # y co-ordinate of lower left corner of middle square\n", " x2 = x + 2/3 * L # x co-ordinate of upper right corner of middle square\n", " y2 = y + 2/3 * L # y co-ordinate of upper right corner of middle square\n", "\n", " # Plot inverted square\n", " plt.fill([x1, x2, x2, x1], [y1, y1, y2, y2], facecolor=\"white\")\n", "\n", " # Repeat for each sub-square\n", " if n > 1:\n", " L = L / 3\n", " sierpinskiCarpet(x, y, L, n - 1) # bottom left square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # bottom square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # bottom right square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # left square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # right square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # top left square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # top square\n", " # sierpinskiCarpet(????, ????, L, n - 1) # top right square\n", "\n", "\n", "# Plot stage 0 square\n", "initialLength = 1\n", "fig, ax = plt.subplots()\n", "plt.fill([0, initialLength, initialLength, 0], [0, 0, initialLength, initialLength], facecolor=\"blue\")\n", "plt.axis(\"equal\")\n", "plt.axis(\"off\")\n", "\n", "# Generate Sierpinsky triangle\n", "numberStages = 2\n", "sierpinskiCarpet(0, 0, initialLength, numberStages)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "UScURMLByFri" }, "source": [ "---\n", "## Fractal Trees\n", "\n", "If you look at trees in nature you may now start to notice that they also have a fractal structure. A main trunk grows upwards until a point where it splits into two or more branches which will also grow for a bit before splitting again. This behaviour continues into the smallest branches begin to sprout leaves. We can attempt to draw a fractal that resembles a tree by doing something similar.\n", "\n", "The Python program below draws a fractal tree up to stage 4. Each branch of this tree splits into two smaller branches, the first branch is 50% of the length of the previous branch and grows at an angle 60 degrees to the left, the second branch is 70% the length of the previous branch and grows at an angle 30 degrees to the right. If we increase the number of stages we get a more detailed tree, although we wouldn't recommend going higher than stage 10 as this will take several minutes to run!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "pSarvyQdyFrj", "scrolled": false }, "outputs": [], "source": [ "def tree(x1, y1, L, angle, n):\n", " # Calculate the end point of the current branch\n", " x2 = x1 + L * np.cos(np.radians(angle))\n", " y2 = y1 + L * np.sin(np.radians(angle))\n", "\n", " if n == 0:\n", " # Add current branch\n", " x.append([x1, x2])\n", " y.append([y1, y2])\n", " if L < 0.3:\n", " colour.append(\"green\")\n", " else:\n", " colour.append(\"brown\")\n", " width.append(10 * L)\n", " else:\n", " # Add branches\n", " tree(x1, y1, L, angle, n - 1) # main branch\n", " tree(x2, y2, 0.5 * L, angle + 60, n - 1) # first branch (50% of the length and 60 degrees to the left)\n", " tree(x2, y2, 0.7 * L, angle - 30, n - 1) # second branch (70% of the length and 30 degrees to the right)\n", "\n", "\n", "# Generate fractal tree\n", "initialAngle = 90\n", "initialLength = 1\n", "numStages = 4\n", "x, y, colour, width = [], [], [], []\n", "tree(0, 0, initialLength, initialAngle, numStages)\n", "\n", "# Plot tree\n", "fig, ax = plt.subplots()\n", "for i in range(len(x)):\n", " plt.plot(x[i], y[i], color=colour[i], linewidth=width[i], solid_capstyle='round')\n", "\n", "ax.axis(\"off\")\n", "ax.set_aspect(\"equal\")" ] }, { "cell_type": "markdown", "metadata": { "id": "dYtjbPWLyFrj" }, "source": [ "---\n", "### Activity 3: Fractal tree\n", "\n", "Amend the Python program below so that it draws a fractal tree where the branches split into three each time with the first branch growing straight up, the second branch growing at an angle of 15 degrees to the left and the third branch 15 degrees to the right and each sub-branch is 75% the length of the previous branch." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "BEXSact-yFrj" }, "outputs": [], "source": [ "def tree(x1, y1, L, angle, n):\n", " # Calculate the end point of the current branch\n", " x2 = x1 + L * np.cos(np.radians(angle))\n", " y2 = y1 + L * np.sin(np.radians(angle))\n", "\n", " if n == 0:\n", " # Add current branch\n", " x.append([x1, x2])\n", " y.append([y1, y2])\n", " if L < 0.3:\n", " colour.append(\"green\")\n", " else:\n", " colour.append(\"brown\")\n", " width.append(10 * L)\n", " else:\n", " # Add branches\n", " tree(x1, y1, L, angle, n - 1) # main branch\n", " # tree(x2, y2, ????, ????, n - 1) # first branch\n", " # tree(x2, y2, ????, ????, n - 1) # second branch\n", " # tree(x2, y2, ????, ????, n - 1) # second branch\n", "\n", "\n", "# Generate fractal tree\n", "initialAngle = 90\n", "initialLength = 1\n", "numStages = 4\n", "x, y, colour, width = [], [], [], []\n", "tree(0, 0, initialLength, initialAngle, numStages)\n", "\n", "# Plot tree\n", "fig, ax = plt.subplots()\n", "for i in range(len(x)):\n", " plt.plot(x[i], y[i], color=colour[i], linewidth=width[i], solid_capstyle='round')\n", "\n", "ax.axis(\"off\")\n", "ax.set_aspect(\"equal\")" ] }, { "cell_type": "markdown", "metadata": { "id": "_o-rYI2hyFrj" }, "source": [ "---\n", "## The Mandlebrot Set\n", "\n", "Perhaps the most famous and beautiful fractal is the [**Mandelbrot set**](https://en.wikipedia.org/wiki/Mandelbrot_set) which was discovered by Polish mathematician Benoit Mandelbrot who was the one who first coined the name *fractal*. The Mandelbrot set is generated by using a very simple formula that is applied to complex numbers, so before we look at how it is generated we are first going to introduce you to complex numbers.\n", "\n", "### Complex numbers\n", "\n", "\n", "\n", "Despite their name, complex numbers are actually quite simple. A complex number is a number that is expressed using two parts: a **real part** and an **imaginary part**. The real part is just a real number, the imaginary part is the **imaginary number** $i = \\sqrt{-1}$ multiplied by a real number. For example, the complex number $z = x + yi$ has a real part $x$ and an imaginary part $y$. We can think of a complex number as being a point on the complex plane where the horizontal axis represents the real part and the vertical axis represents the imaginary part. The distance of the point from the origin is known as the **absolute value** (or modulus) of $z$ and is denoted by $|z|$. Using Pythagoras' theorem we calculate the absolute value of a complex number using\n", "\n", "\\begin{align*}\n", " |z| = \\sqrt{x^2 + y^2}.\n", "\\end{align*}\n", "\n", "To add two complex numbers we simply add the real parts and the imaginary parts separately, for example\n", "\n", "\\begin{align*}\n", " (x + yi) + (a + bi) = (x + a) + (y + b)i.\n", "\\end{align*}\n", "\n", "To multiply of two complex numbers we multiply out two bracketed expressions using the FOIL technique (first, outside, inside last), for example\n", "\n", "\\begin{align*}\n", " (x + yi) (a + bi) &= xa + xbi + yai + ybi^2 = (xa - yb) + (xb + ya)i.\n", "\\end{align*}\n", "\n", "Since squaring a number is multiplying a number by itself the square of a complex number is\n", "\n", "\\begin{align*}\n", " (x + yi)^2 = (x + yi)(x + yi) = x^2 - y^2 + 2xyi.\n", "\\end{align*}\n", "\n", "\n", "### Generating the Mandelbrot set\n", "\n", "The Mandelbrot set is generated by the simple equation\n", "\n", "$$z_{n+1} = z_n^2 + c,$$\n", "\n", "where $z_n$ is the current number in a sequence of numbers, $z_{n+1}$ is the next number, $z_0=0$ is the starting number and $c = x + yi$ is some complex number. The first time we calculate this equation we have $z_1 = c$, the second time we have $z_2 = c^2 + c$, the third time we have $z_3 = (c^2 + c)^2 + c$ and so on. This method of repeating a calculation the same equation is known as **iteration**. We continue to iterate this equation until either the absolute value of $z_n$ is greater than 2 or a maximum number of iterations has been reached. If $|z_n|>2$ then we say that the number $c$ has *escaped* and is not a member of the Mandelbrot set, else if at the end of the iterations the absolute value of $z_n$ is still less than 2 we say the number $c$ is a member of the Mandelbrot set and we mark this point on the complex plane. If we repeat this for other points in the complex plane a beautifully complex shape starts to appear.\n", "\n", "For example consider the number $c = 1$, the iterations are:\n", "\n", "\\begin{align*}\n", " z_1 &= 1, & |z_1| &= 1, \\\\\n", " z_2 &= 1^2 + 1 = 2, & |z_2| &= 2, \\\\\n", " z_3 &= 2^2 + 1 = 5, & |z_3| &= 5.\n", "\\end{align*}\n", "\n", "So after 3 iterations the number $c=1$ has escaped so is not a member of the Mandelbrot set. Let's try the number $c = i$:\n", "\n", "\\begin{align*}\n", " z_1 &= i, & |z_1| &= 1, \\\\\n", " z_2 &= i^2 + i = -1 + i, & |z_2| &\\approx 1.41, \\\\\n", " z_3 &= (-1 + i)^2 + i = -i, & |z_3| &= 1, \\\\\n", " z_4 &= (-i)^2 + i = -1 + i, & |z_4| &\\approx 1.41.\n", "\\end{align*}\n", "\n", "Here $z_4=z_2$ so the values of $z_n$ will cycle between $-1+i$ and $-i$. This means $|z_i|$ will always be less than 2 so $c = i$ is a member of the Mandelbrot set.\n", "\n", "The Python program below draws the Mandelbrot set for a region on the complex plane defined by `x1` and `y1` which are the co-ordinates of the bottom-left corner of the region and `width` and `height` which define the size of the region. `numXPixels` is the number of points in the $x$-direction and `maxIterations` which is the maximum number of iterations that are calculated. Run the code cell to see the result." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "LwfheJqZyFrk" }, "outputs": [], "source": [ "def mandelbrot(x1, y1, width, height, numXPixels, maxIterations):\n", "\n", " # Define C, Z and M arrays\n", " numYPixels = int(numXPixels * height / width)\n", " X, Y = np.meshgrid(np.linspace(x1, x1 + width, numXPixels), np.linspace(y1 + height, y1, numYPixels))\n", " C = X + Y * 1j\n", " Z = np.zeros(C.shape).astype(complex)\n", " M = np.zeros(C.shape)\n", "\n", " # Perform iterations\n", " for n in range(maxIterations):\n", " i = np.where(abs(Z) < 2)\n", " Z[i] = Z[i] ** 2 + C[i]\n", " M[i] += 1\n", "\n", " # Find the points that are still inside the Mandelbrot set\n", " M = M == maxIterations\n", "\n", " # Plot the Mandelbrot set\n", " fig, ax = plt.subplots(figsize=(8, 8 * height / width))\n", " plt.imshow(M, extent=[x1, x1 + width, y1, y1 + height], cmap=cmap)\n", " plt.xlabel(\"real(c)\")\n", " plt.ylabel(\"imag(c)\")\n", "\n", "\n", "# Define Mandelbrots set parameters\n", "x1, y1, width, height, numXPixels, maxIterations = -2.25, -1.25, 3, 2.5, 400, 20\n", "\n", "# Plot the Mandelbrot set in black and white\n", "cmap = colors.ListedColormap(['white', 'black'])\n", "mandelbrot(x1, y1, width, height, numXPixels, maxIterations)" ] }, { "cell_type": "markdown", "metadata": { "id": "ZmOExcl4yFrk" }, "source": [ "### Adding colour to the Mandelbrot set\n", "\n", "Some numbers which are not members of the Mandelbrot set take longer to escape than others. We can represent this on our Mandelbrot set by plotting each point using a colour which depends on its escape count. Run the program below to plot the Mandelbrot set in colour." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Alz1yaGUyFrk" }, "outputs": [], "source": [ "import matplotlib\n", "\n", "def mandelbrot(x1, y1, width, height, numXPixels, maxIterations):\n", "\n", " # Define C, Z and M arrays\n", " numYPixels = int(numXPixels * height / width)\n", " X, Y = np.meshgrid(np.linspace(x1, x1 + width, numXPixels), np.linspace(y1 + height, y1, numYPixels))\n", " C = X + Y * 1j\n", " Z = np.zeros(C.shape).astype(complex)\n", " M = np.zeros(C.shape)\n", "\n", " # Perform iterations\n", " for n in range(maxIterations):\n", " i = np.where(abs(Z) < 2)\n", " Z[i] = Z[i] ** 2 + C[i]\n", " M[i] += 1\n", "\n", " # Smooth out the escape count\n", " i = np.where(M < maxIterations)\n", " M[i] += 1 - np.log(np.log2(abs(Z[i])))\n", "\n", " # Plot the Mandelbrot set\n", " fig, ax = plt.subplots(figsize=(8, 8 * height / width))\n", " plt.imshow(M, extent=[x1, x1 + width, y1, y1 + height], cmap=cmap)\n", " plt.xlabel(\"real(c)\")\n", " plt.ylabel(\"imag(c)\")\n", "\n", "\n", "# Define Mandelbrots set parameters\n", "x1, y1, width, height, numXPixels, maxIterations = -2.25, -1.25, 3, 2.5, 400, 40\n", "\n", "# Define colormap\n", "jet = matplotlib.colormaps.get_cmap(\"jet\") # get the colormap\n", "colours = jet(np.linspace(0, 1, 256)) # get the colours from the colormap\n", "colours = np.vstack((colours, [0, 0, 0, 1])) # add a black layer to the colormap\n", "cmap = colors.ListedColormap(colours) # save colormap\n", "\n", "# Plot the Mandelbrot set\n", "mandelbrot(x1, y1, width, height, numXPixels, maxIterations)" ] }, { "cell_type": "markdown", "metadata": { "id": "pbB7l42IyFrk" }, "source": [ "---\n", "### Activity 4: Exploring the Mandelbrot set\n", "\n", "The amazing thing about the Mandelbrot set is that when we zoom in to areas surrounding the points in the Mandelbrot set we see that these form complicated and beautiful shapes. There are even small Mandelbrot sets found which we can also zoom into so the Mandelbrot set has the self-similarity property found in other fractals. The code in the code cell below plots the region of the Mandelbrot set defined by bottom-left and top-right co-ordinates (-1.05, 0.25) and (-0.95, 0.35) and calculated for 200 iterations. Explore the Mandelbrot set by changing these values (increase `maxIterations` to get more detail)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Jr4BNkQlyFrk", "scrolled": false }, "outputs": [], "source": [ "# Plot the mandelbrot set (inputs: x1, y1, width, height, numXPixels, maxIterations)\n", "mandelbrot(-1.05, 0.25, 0.1, 0.07, 400, 300)" ] }, { "cell_type": "markdown", "metadata": { "id": "PNQbml8PyFrl" }, "source": [ "### Zooming into the Mandlebrot set\n", "\n", "The video below was created by YouTube creator [Maths Town](https://www.youtube.com/channel/UC6qEdtxp_IAaVrNAHUIhHbQ) and zooms into the Mandelbrot set showing the infinite complexity of the fractal." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "L5lyCDf0yFrl", "scrolled": true }, "outputs": [], "source": [ "from IPython.display import YouTubeVideo\n", "YouTubeVideo('pCpLWbHVNhk', width=600)" ] }, { "cell_type": "markdown", "metadata": { "id": "-547QpJRyFrl" }, "source": [ "---\n", "## Useful Links\n", "\n", "If you would like to further explore Python and fractals you may find the following links useful:\n", "\n", "- [Anaconda](https://www.anaconda.com/products/distribution) - a suite of software tools that includes Jupyter Notebook and Python. Download and install on your computer to write and run Jupyter notebooks\n", "- [Google Colab](https://colab.research.google.com/) - run Jupyter notebooks in the cloud using Google Colab (you will need to have a Google account to do this)\n", "- [Fractals](https://en.wikipedia.org/wiki/Fractal)\n", "- [Koch snowflake](https://en.wikipedia.org/wiki/Koch_snowflake)\n", "- [Sierpinski Triangle](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_triangle)\n", "- [Sierpinski carpet](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_carpet)\n", "- [Fractal tree](https://en.wikipedia.org/wiki/Fractal_canopy)\n", "- [Mandelbrot set](https://en.wikipedia.org/wiki/Mandelbrot_set)\n", "\n", "© Dr Jon Shiach 2023" ] } ], "metadata": { "colab": { "provenance": [], "include_colab_link": true }, "interpreter": { "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.3" } }, "nbformat": 4, "nbformat_minor": 0 }